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ABSTRACT 


The  performance  of  the  basic  self- tuning  regulator  (STR)  for 
control  of  top  product  composition  of  a  binary  distillation 
column  is  studied  experimentally  and  by  simulation.  The 
self-tuning  regulator  provides  satisfactory  control  behavior 
for  step  changes  in  feed  flow  rate  and  set  point  even  though 
the  regulator  is  not  formulated  for  such  changes.  A  PI 
controller,  once  well -tuned  provides  excellent  performance 
and  outperforms  the  STR,  when  used  with  zero  initial 
parameter  estimates. 

In  the  case  of  bottom  product  composition  control, 
simulations  have  shown  that  due  to  the  long  time  delay  in 
the  bottom  composition  analysis,  the  PI  controller  is  not 
able  to  provide  satisfactory  control.  The  STR,  with  the 
advantage  of  an  inherent  predictor  in  the  model  formulation, 
manages  to  achieve  a  much  improved  control  performance  once 
reasonable  parameter  estimates  were  obtained  with  the 
multiple  pass  method.  In  this  work,  the  principle  STR 
design  parameters  ^>o  and  P(0)  have  been  shown  to 
significantly  effect  the  performance  of  the  regulator. 
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1.  INTRODUCTION 


1.1  Introduct ion 

The  distillation  column,  despite  its  extensive  use  in 
industry,  remains  among  the  most  difficult  unit  operations 
to  control.  The  physical  complexity  of  the  distillation 
process  and  the  inherent  system  nonlinearity  combine  to 
provide  a  major  challenge  for  those  responsible  for 
designing  control  systems. 

With  the  availability  of  process  control  computers, 
different  advanced  control  schemes  have  constantly  been 
evaluated  as  part  of  a  continuing  research  activity  in  hope 
of  achieving  better  efficiency,  higher  throughput,  lower 
operating  cost  or  any  attractive  combination  of  these.  The 
conventional  three  term  controller,  even  properly  tuned  for 
the  system's  dynamic  behavior  is  supposedly  a  compromise 
between  disturbance  minimization  and  set  point  changes.  It 
has  often  been  criticized  for  not  yielding  an  optimal 
control  and  for  its  inability  to  handle  nonlinear  dynamic 
behavior . 

To  design  a  suitable  controller,  an  accurate  model  is 
usually  required.  Among  advanced  control  techniques,  those 
based  on  linear  state  space  models  have  not  been  very 
successful  due  to  the  difficult  problem  of  representing  the 
dynamic  behavior  of  the  process  system  using  a  reduced  state 
space  model  (1).  For  feedforward  control,  even  if  a  good 
model  is  obtained,  the  process  operating  conditions  may 
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change  and  the  model  becomes  inadequate  (2,3). 

An  adaptive  control  strategy,  developed  by  Astrom  and 
WittenmarK  appears  to  overcome  these  difficulties  .  It 
allows  on-line  tuning  of  the  model  parameters.  This 
strategy  has  been  termed  as  the  self-tuning  regulator  (4,5). 

The  basic  self-tuning  regulator  is  designed  for  the 
regulator  problem.  It  is  designed  to  control  single  input, 
single  output  systems  with  unknown  but  constant  parameters. 
The  regulator  can  also  be  applied  to  systems  with  slowly 
varying  parameters.  The  algorithm  combines  a  recursive 
least  squares  estimator  with  a  minimum  variance  regulator- 
computed  from  the  estimated  model. 


1.2  Objective  of  Thesis: 

This  work  is  directed  toward  assessing  the 
effectiveness  of  the  basic  self-tuning  regulator  for  control 
of  a  binary  distillation  column.  By  applying  the  technique, 
a  better  understanding  of  its  relative  advantages, 
di sadvantages  and  applicability  in  the  process  industry  will 
be  achieved. 

In  this  investigation,  the  performance  of  the 
self- tuning  regulator  for  control  of  top  product  composition 
is  studied  experimentally  and  by  simulation.  The  control  of 
bottom  composition  has  been  studied  by  simulation  for  load 
and  set  point  changes.  The  resulting  control  performance  is 
compared  with  that  of  a  well-tuned  conventional  proportional 
plus  integral  (PI)  controller. 
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1.3  Outline  of  Thesis 

The  thesis  is  organized  in  four  parts: 

The  first  section,  consisting  of  Chapter  2,  presents 
the  mathematical  formulation  of  the  self-tuning  regulator. 

A  literature  survey  is  not  included  since  detail  surveys 
have  been  given  by  Chang  (6)  and  Lieuson  (7). 

Chapter  3,  which  constitutes  the  second  section, 
provides  the  necessary  process  information  for  the 
application  of  the  adaptive  control  schemes. 

The  third  section,  Chapters  4  and  5  give  the  results  of 
control  of  top  composition  and  bottom  composition 
respectively. 

The  last  section,  Chapter  6  summarizes  the  overall 
conclusions  concerning  the  utilization  of  the  selftuning 
regulator . 

Note:  The  work  of  this  thesis  was  completed  during 
1977.  However,  the  thesis  was  submitted  for  examination  in 
September  1980. 
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2.  THE  SELF- TUNING  REGULATOR 


2.1  Theoretical  Formulation 

The  self-tuning  regulator  of  Astrom  and  Wittenmark  (4) 
was  derived  using  the  following  single- input ,  single-output 
linear  stochastic  model  of  the  process  as  a  starting  point: 

y(t)  +  aly(t-1)  +  ...  +  any(t-n)  =  bt  u ( t -K- 1 )  +  ...  + 
bnu(t-k-n)  +Xle(t)  +  cAe(t-1)  +  ...  +  c^e(t-n)]  (2.1) 

where : 

y  is  the  normalized  output  deviation  variable 
u  is  the  normalized  input  deviation  variable 
t  denotes  the  sampling  instant 

(e(t),t  =  0,1,2,...)  is  a  sequence  of  independent  normal 
(0,1)  random  variables  which  represent  the  net 
effect  of  external  disturbances, 
k  is  a  nonnegative  integer  which  expresses  the  process 
time  delay  as  a  multiple  of  the  sampling  interval. 

{ a,.  }  f  {b(  }  ,  { c  i. }  and  X  are  time  invariant 
parameters . 

The  stochastic  model  in  Equation  (2.1)  has  received 
considerable  attention  not  only  in  the  stochastic  control 
theory  literature  but  also  by  workers  in  time  series 
analysis  (8,9).  This  model  can  also  be  expressed  in 
backward  shift  operator  notation  as: 

A(q”*  )y(t)  =  B  (q**1  )u(  t-k- 1  )  +  C(q  1  )e(t)  (2.2) 
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where  q”*  denotes  the  backward  shift  operator;  i.e.  q’*x(t) 
=  x ( t ■ 1 )  and  polynomials  A,  B  and  C  are  defined  by: 

A ( q” 1  )  =  1  +  atq-1  +  a^q  2  +  ...  +  anq'n  (2.3) 

-1  -l  -rv+*l 

B  (q  )  =  bt  +  b2q  +  ...  +  bnq  (2.4) 

C(q'1)  =  1  +  c^q'1  +  ...  +  cnq  n  (2.5) 

The  minimum  variance  regulator  minimizes  the  loss 
function,  V  =  E{y  (t)}  where  E  denotes  the  expectation 
operator.  Astrom  (4)  has  derived  the  following  minimum 


variance  control  law: 

u(t)=  -G(ql)/B(q”1)F(q1)*y(t)  (2.6) 

where  the  polynomials 

Flq'1  )  =  1  +  fjq'1  +  .  .  .  +  fkqk  (2.7) 

G(q_1  )  =  ga  +  ...  +  g„_1q'n'vl  (2.8) 

are  related  to  A(q’*  )  and  C(q’*  )  by: 

C(q_1)  =  A(q'1)F(q"1)  +  q~k-1  G(q"1)  (2.9) 


Since  the  model  parameters  {a*},  { b ^ } ,  {c^l  and  \  are 
unknown  and  constant,  the  self-tuning  strategy  consists  of 
estimating  the  model  parameters  on-line  and  basing  the 
control  calculation  at  each  sampling  instant  on  the 
estimated  parameters.  A  wide  variety  of  SIR  algorithms  can 
be  derived  depending  on  the  particular  estimation  algorithm 
and  control  strategy  that  are  selected  (5).  In  this 
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investigation,  a  recursive  linear  least  squares  estimator  is 
used  in  conjunction  with  a  minimum  variance  controller  as 
originally  suggested  by  Astrom  and  Wittenmark  (4). 

In  principle,  the  least  squares  calculations  could  be 
based  on  updating  the  parameters  of  the  model  in  Equation 
(2.1)  but  this  would  involve  a  complicated  estimation 
procedure  which  is  not  suitable  for  on-line  computation.  By 
assuming  C;.  =  0  for  i  =  1,2,...n  the  parameter  estimation 
problem  reduces  to  a  conventional  linear  least  squares 
problem  which  can  be  easily  formulated  as  a  recursive 
algorithm.  The  on-line  control  calculations  can  be 
simplified  by  i ncorporat i ng  the  identity  in  Equation  (2.9) 
into  the  model  structure,  resulting  in  the  following 
"predictive"  model: 

y(t)  +  P'4  y  ( t -k- 1  )  +  ...  +P^y(t-k-n)  = 

j^>  [u(  t-k-  1  )  +  ...+  u(  t-L-K-  1 )  ]  +  £(t)  (2.10) 

where  L  =  n  +  K  -  1  and  the  disturbance  fc. ( t )  is  a  moving 
average  of  order  k  of  the  driving  noise,  e(t). 

If  the  coefficients  in  the  predictive  model  are  known 
and  constant,  then  the  minimum  variance  regulator  for  the 
predictive  model  in  Equation  (2.10)  is: 

u ( t )  =  1/j^0[^y(t)  +  .  .  .•  +  oCn -yCt-n+1)] 

-^u(t-l)  -  ...  -£Lu(t-L)  (2.11) 

In  the  self-tuning  regulator,  {0O}  and  { ^\  }  are 

A  A 

estimated  on-line  and  the  estimates  {c>6  }  and  {3*}  are  used 
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in  the  minimum  variance  controller  of  Equation  (2.11). 
However,  it  may  not  be  possible  to  accurately  estimate  all 
of  the  parameters  during  closed-loop  operation  since  the 
inputs  are  correlated  with  the  disturbance  {£(£')}•  To  avoid 
this  difficulty,  Astrom  and  Wittenmark  (4)  suggested  that  ^>o 
should  be  held  at  some  specified  value.  For  constant  , the 
STR  control  algorithm  is: 

u(t)  =  1/|>0  ( t  )y(  t )  +  ...  +  c/n  y(  t-n+1 )  ] 

-  (t)u(t-l  )  "...  -^.(t)u(t-L)  (2.12) 

which  can  be  written  as: 

u  ( t )  =  -1A  e(t)J(t)  (2.13) 

where 

^T(t)  =  [  -y(t),  -  y ( t - 1 ) ,  -  ...  -  y(t-n+1), 

+  J^O  u(t-1  )  ,  +  ...  ,+  ^,u(t-L)  ]  (2.14) 

and 

e(t)  =  [ ofj  ( t ) . Uz(t),<Z(t),  ^(t) . |L(t)i  (2.15) 

A 

The  parameter  estimates,  0(t)  are  calculated  using 
recursive  least  squares  identification  with  exponential 
weighting  of  past  data  (10).  The  algorithm  is  given  by: 

0(t)  =  0 ( t- 1 )  + 

K(t)  (y(  t )  -0ou(t-k-1)  -ir(t-K-1)|(t-1))  (2.16) 

4 

K  ( t )  =  P(t-1)+(t-k-1)[1  +fT(t-k-1)p(t-1)t(t-k-1)r  (2.17) 
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P(t)  =  1/u  {  P ( t-  1 ) 

-  K(t)[  1  +  ^T(  t -K-  1  )  P  ( t  - 1  )^(  t -K- 1 )  k7(  t )  J  }  (2.18) 

is  the  estimated  value  of  6(t) 
is  the  gain  vector  {  (n+1)  x  1  } 
is  proportional  to  the  covariance 
matrix  {  (n+1  )  x  (n+1 ) } 
is  the  exponential  forgetting  factor 

In  the  estimatation  phase,  K(t)  is  calculated  first, 
followed  by  £(t)  and  0(t)  as  discussed  in  Appendix  A. 


where 

0(t) 

K(t) 

P(t) 


2.2  Parameter  Selection 

Before  implementing  the  self-tuning  regulator,  the 
constants  and  initial  values  listed  below  must  be  specified: 
model  order  n 
time  delay  K 
scaling  factor  j?>0 
exponential  forgetting  factor  fx 
initial  parameter  estimates  0(0) 
initial  covariance  matrix  £(0) 
constraint  limits  on  control  output  U.L 
Astrom  and  WittenmarK  (4)  and  Chang  (6)  have  presented 
guidelines  for  selecting  appropriate  values  for  the  above 
constants.  These  recommendations  have  been  used  in  this 


work . 
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2.2.1  Model  Order,  n: 

In  order  that  the  process  be  adequately  represented, 
the  self- tuning  regulator  must  have  enough  parameters  so  it 
could  converge  to  the  optimal  minimum  variance  regulator. 

For  the  zero  initial  estimate  case,  if  a  higher  order  model 
is  used,  it  will  need  more  sampling  intervals  before  the 
information  vector  (  +  (t))  and  the  parameter  vector  (0(t)) 
have  non-zero  elements.  Consequently,  this  will  result  in 
poorer  control  during  the  initial  transient  period. 

2.2.2  System  Time  Delay,  k: 

The  choice  of  k,  the  nonnegative  integer  to  represent 
the  time  delay  of  the  process  in  terms  of  an  integral  number 
of  sampling  intervals  is  extremely  important-.  It  directly 
affects  the  number  of  parameters  to  be  estimated. 
Overestimating  this  parameter  will  result  in  a  sluggish 
control  performance  whereas  underestimating  k  will  cause 
cycling  or  bang  bang  control. 

2.2.3  Scaling  Factor,  : 

The  value  selected  for  should  be  close  to  the 
actual  value  for  the  system  to  be  controlled.  From  the 
consideration  of  stable  operation,  if  the  value  used  is 
much  lower  than  the  true  value,  it  may  cause  the  parameters 
to  diverge  (4).  The  effect  of  the  choice  of  the  on  the 
controller  performance  is  investigated. 
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2.2.4  Exponential  Forgetting  Factor, jx  : 

The  use  of  a  weighting  factor  less  than  1  allows  the 
weighting  out  of  previous  data.  This  feature  enables 
tracking  of  slow  time  varying  parameters.  However,  in  this 
study,  due  to  the  short  duration  of  the  runs,  the  parameters 
are  not  considered  to  be  time  varying,  so  as  a  result,  j\  is 
set  to  1. 

2.2.5  Initial  Parameter  Estimates,  0(0): 

The  initial  estimates  will  frequently  be  set  to  zero 
when  there  is  no  a  priori  information.  The  use  of  other 
than  zero  initial  values  for  the  second  pass  or  subsequent 
runs  is  investigated  in  this  study. 

2.2.6  Initial  Covariance  Matrix,  £(0): 

The  initial  covariance  matrix  used  reflects  the 
confidence  in  the  initial  parameter  estimates.  A  higher 
£(0)  value  should  increase  the  convergence  rate  of  the 
parameters  at  each  sampling  interval.  The  effect  on 
controller  performance  of  this  parameter  is  investigated  in 
this  study. 

2.2.7  Limits  on  control  signal  ampl i tude,  UL  : 

Control  signal  limits  are  implemented  to  safeguard 
against  run  away  conditions.  However,  too  tight  a  limit  may 
slow  the  convergence  rate  of  the  parameters. 
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2.3  APPLICATION  OF  THE  STR  ALGORITHM 

In  the  past  years,  some  of  the  industrial  processes 
that  have  been  successfully  controlled  by  the  self-tuning 
regulator  include  a  paper  machine  (11),  an  ore  crusher  (12), 
a  heat  exchanger  (13)  and  a  super  tanker  (14).  At  the 
University  of  Alberta,  Chang  (6)  has  studied  the  self-tuning 
regulator  control  of  a  double  effect  evaporator  and  Sastry 
et  a  1 . ( 1 5 )  investigated  top  product  composition  control  on  a 
distillation  column. 

For  implementing  the  self-tuning  regulator,  limited  a 
priori  knowledge  of  the  process  dynamics  and  noise 
character i st i cs  of  the  system  to  be  controlled  are  needed. 

The  moderate  computation  requirement  of  the  algorithm 
makes  it  very  attractive  for  industrial  use.  As  the 
algorithm  is  adaptive,  it  is  ideal  for  use  in  dealing  with 
processes  that  have  unknown  process  dynamics,  nonlinearity, 
or  changes  in  operating  conditions.  The  regulator  is  best 
suited  for  a  situation  where  the  control  objective  is  to 
minimize  output  deviation. 

A  possible  drawback,  however,  would  be  the  regulator's 
extreme  sensitivity  to  parameter  variation  for  nonminimum 
phase  systems  (4).  The  necessity  of  extending  the  regulator 
to  include  feedforward  compensation  for  measurable 
disturbances  and  set  point  following  was  recognized  by 
Wittenmark  (16).  Clarke  and  Gawthrop  (17)  have  recommended 
the  use  of  a  cost  function  that  incorporates  system  input, 
output  and  set  point.  In  Clarke's  algorithm,  control  output 


, 


is  reduced  at  the  expense  of  extra  output  variance. 

In  this  study,  the  basic  self-tuning  regulator,  as 
proposed  by  Astrom  is  evaluated  on  a  pilot  scale  binary 
distillation  col umn . 
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3.  PROCESS  CHARACTERIZATION 


3.1  Process  Equipment 

The  experimental  studies  were  performed  on  a  pilot 
scale  distillation  column  interfaced  to  an  IBM  1800  data 
acquisition  and  control  computer.  A  schematic  diagram  of 
the  equipment  is  shown  in  Figure  3.1.  The  22.86  cm. 
diameter  column,  containing  eight  bubble  cap  trays  (four 
caps  per  tray)  has  a  30.48  cm.  tray  spacing.  The  column  is 
equipped  with  a  total  condenser  and  thermosyphon  type 
reboiler.  The  column  operates  at  atmospheric  pressure,  with 
the  feed  at  48%  by  weight  methanol  entering  the  fourth  tray. 
Top  product  composition  is  96.30%  by  weight  methanol  and  the 
bottom  product  composition  is  5.50%  by  weight  methanol.  The 
top  composition  is  measured  by  means  of  an  in-line 
capacitance  probe.  The  bottom  composition  is  measured  with 
a  Hewlett-Packard  model  5720A  gas  chromatograph  using  an 
in-line  liquid  sampling  system.  The  chroma togr arris  are 
analyzed  on  a  two  hundred  fifty-six  second  cycle.  All 
measurements  are  transmitted  to  the  IBM  1800  for  data 
logging  and/or  control.  Further  specific  details  concerning 
the  column  are  given  in  the  theses  by  Svrcek  (18)  and  Pacey 
(19), 

The  experimental  operating  conditions  used  in  this 
work,  which  differ  somewhat  from  those  of  previous  work,  are 
given  in  Table  3.1. 
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Figure  3.1  Schematic  Diagram  of  the  Distillation  Column 
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Table  3.1  Experimental  Operating  Conditions  (Run  El) 


Steam  Flow 


Reflux  Flow 


Feed  Flow 


18.0  g/s  Feed  Composition  48.0  wt% 
10.3  g/s  Top  Product  Comp.  96.3  wt% 
11.5  g/s  Bottom  Product  Comp.  5.5  wt% 


Top  Product  Flow  8.4  g/s  Feed  Temperature  74.1  C 

Bottom  Product  Flow  9.6  g/s  Reflux  Temperature  64.2  C 


3.2  Process  Model 

A  nonlinear  ordinary  differential  equation  model,  based 
on  the  mass  and  energy  balance  relationships  for  each  of  the 
trays,  reboiler  and  condenser  is  utilized  for  simulating  the 
dynamic  behavior  of  the  column  (20).  The  initial  analysis 
results  in  30  differential  equations  containing  140 
variables  but  the  number  of  unknown  variables  can  be  reduced 
to  20  by  specifying  certain  inputs  and  outputs  and  by  making 
the  following  assumptions: 

1.  Constant  heat  loss  from  each  stage 

2.  Constant  tray  efficiency 

3.  Constant  liquid  holdup 

4.  Negligible  vapor  holdup 

Details  of  the  model  development  and  numerical  solution 
are  contained  in  the  thesis  of  Simonsmeier  (20).  The  steady 
state  operating  conditions  used  for  simulation  ,  which  were 
similar  to  those  used  experimentally  are  given  in  Table  3.2 


- 
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Table  3.2  Operating  Conditions  For  Simulation 


Top  Product  Flow 


Steam  Flow 


Reflux  Flow 


Feed  Flow 


18.0  g/s  Feed  Composition  48.0  wt% 
10.3  g/s  Top  Product  Comp.  95.2  wt% 
11.75  g/s  Bottom  Product  Comp.  5.5  wt% 
8.6  g/s  Feed  Temperature  74  C 


Bottom  Product  Flow  9.4  g/s  Reflux  Temperature  64  C 

3.2.1  Dynamic  Testing 

The  objective  of  the  experimental  work  was  the  control 
of  either  only  the  top  composition  or  bottom  composition. 
Control  was  evaluated  for  load  and  set  point  changes,  i.e. 
for  both  regulatory  and  servo  control . 

The  main  disturbances  to  a  column  are  usually  changes 
in  feed  flow  rate  and  the  feed  composition.  The  feed  flow 
rate  tends  to  vary  much  faster  and  more  often  than  does  the 
feed  composition  and  in  fact,  both  Shinskey  (21)  and  Lupfer 
(22)  suggest  that  feedback  control  designed  on  the  basis  of 
feed  flow  rate  disturbances  should  be  sufficient  to 
compensate  for  the  gradual  feed  composition  changes.  Hence 
feed  flow  rate  was  chosen  as  the  system  load  disturbance. 

Steam  flow  has  normally  been  reserved  for  control  of 
bottom  composition.  By  following  typical  industrial 
practice,  reflux  flow  rate  rather  than  top  product  flow  rate 
was  chosen  as  the  manipulated  variable  for  top  product 
composition  control.  Top  product  flow  was  used  to  control 
the  condenser  level. 
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In  order  to  implement  the  self-tuning  regulator 
algorithm,  it  is  necessary  that  a  predictive  model  be 
specified.  Time  delay  in  the  process  and  the  number  of 
parameters  in  the  regulator  thus  have  to  be  determined  in 
advance.  In  addition,  high  and  low  limits  for  the  control 
signal  must  be  chosen.  To  provide  reliable  information  for 
choice  of  these  values,  dynamic  testing  of  the  column  by 
using  open  loop  step  disturbance  runs  were  performed.  Feed 
flow  rate  step  changes  of  ±  20%,  from  the  steady  state 
value,  were  introduced  to  study  the  response  of  the  column, 
while  the  two  manipulated  variables,  steam  and  reflux  flows 
were  changed  ±  5%  and  ±  10%  respectively.  The  responses  for 
the  step  changes  are  plotted  in  Figures  3.2  to  3.7.  The 
arrow  on  the  abcissa  scale  in  each  plot  indicates  the  time 
at  which  the  disturbance  is  introduced. 

Examination  of  the  responses  plotted  in  Figures  3.2  to 
Figure  3.7  shows  that  reasonable  character izat ion  of  the 
different  open  loop  responses  is  possible  using  a  first 
order  plus  time  delay  transfer  function  although  the  values 
of  the  gain  and  time  constant  are  found  to  be  strongly 
dependent  upon  the  direction  and  magnitude  of  the 
disturbances.  The  initial  and  final  values  of  top  and 
bottom  composition  for  the  different  step  changes  in  the 
disturbances  are  summarized  in  Table  3.3 
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Figure  3.2  Dynamic  Response  of  the  Terminal  Compositions  to 

a  Feed  Flow  Rate  Step  Change  of  +20% 
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Figure  3.3  Dynamic  Response  of  the  Terminal  Compositions  to 

a  Feed  Flow  Rate  Step  Change  of  -20% 
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Figure  3.4  Dynamic  Response  of  the  Terminal  Compositions  to 

a  Reflux  Flow  Rate  Step  Change  of  +10% 
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Figure  3.5  Dynamic  Response  of  the  Terminal  Compositions  to 

a  Reflux  Flow  Rate  Step  Change  of  -10% 
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Figure  3.6  Dynamic  Response  of  the  Terminal  Compositions  to 

a  Steam  Flow  Rate  Step  Change  of  +5% 
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Figure  3.7  Dynamic  Response  of  the  Terminal  Compositions  to 
a  Steam  Flow  Rate  Step  Change  of  -5% 
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Table  3.3  Summary  Of  Product  Composition  Changes 

For  Step  Changes  In  Feed,  Reflux  and 
Steam  Flow  Rates 


Run 

Type  Of 

Di s turbance 

Initial 

XD 

Values 

XB 

Final 

XD 

Va 1 ues 

XB 

E2 

+20%  Feed 

96.4 

4.5 

96.4 

17.8 

E3 

-20%  Feed 

96.4 

4.5 

93.6 

0.3 

E4 

+10%  Reflux 

96.2 

5.2 

97.0 

9.5 

E5 

-10%  Reflux 

96.3 

4.7 

95.3 

1 .  1 

E6 

+5%  Steam 

96.4 

4.5 

95.3 

1  .  1 

E7 

-5%  Steam 

96.2 

4.8 

96.7 

10.8 

3.2.2  Nonlinear  Model 

The  simulated  responses  of  the  top  and  bottom 
composition  obtained  from  the  nonlinear  model  for  feed  flow 
rate  step  changes  of  ±  20%  are  compared  to  the  experimental 
column  responses  in  Figures  3.8  and  3.9.  The  output 
deviation  responses  from  the  nonlinear  model  agree 
reasonably  closely  with  the  experimental  column  outputs. 

In  examining  the  column  output  responses,  it  can  be 
seen  that  the  system  exhibits  highly  nonlinear  behavior. 

For  the  20%  decrese  in  feed  flow  rate,  the  top  composition 
decreases  by  approximately  three  weight  percent  methanol , 
whereas  for  the  same  magnitude  feed  change  in  the  opposite 
direction  there  is  no  effect  on  the  top  composition  at  all. 
Since  the  top  composition  proved  to  be  insensitive  to 
increases  in  feed  flow  rate,  no  control  runs  were  performed 
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Figure  3.8  Comparison  of  Predicted  and  Experimental  Dynamic 

Response  of  the  Terminal  Compositions  to  a  Feed 
Flow  Rate  Step  Change  of  +  20% 
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Figure  3.9  Comparison  of  Predicted  and  Experimental  Dynamic 
Response  of  the  Terminal  Compos i t ions  to  a  Feed 
Flow  Rate  Step  Change  of  -20% 
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for  the  step  increase  in  feed  flow  rate. 

3.2.3  Transfer  Function  Model 

As  mentioned  in  section  3.2.1,  a  first  order  plus  time 
delay  transfer  function  provides  a  reasonably  good 
representation  of  column  dynamic  behavior  for  the  conditions 
used.  This  is  important  since  the  choice  of  K,  which 
represents  the  time  delay  in  the  STR  model,  is  crucial.  Too 
large  a  value  of  K  can  give  an  output  variance  which  is 
larger  than  necessary.  Also,  it  can  prove  to  be  difficult 
to  get  the  estimates  to  converge  at  times.  Too  small  a  K 
makes  it  very  difficult  to  get  good  control  and  may  even 
make  it  impossible  to  stabilize  the  system  (16). 

3.2.3. 1  Top  Composition  Time  Delay 

Analysis  of  the  experimental  responses  in  section 
3.2.1,  indicates  that  the  time  delay  for  the  top  product 
composition,  (XD  vs  Reflux)  is  nearly  constant  at  one 
minute.  Therefore  with  a  control  cycle  of  sixty-four 
seconds,  k  was  chosen  to  be  one. 

3. 2. 3. 2  Bottom  Composition  Time  Delay 

The  dynamic  responses  of  the  bottom  composition  to 
steam  flow  rate  step  changes  of  ±  5%  provide  a  good 
indication  of  the  time  delay,  (XB  vs  Steam).  The  total  time 
delay  is  estimated  to  be  between  five  to  six  sampling 
intervals.  The  system  time  delay  is  due  to 
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(i)  The  gas  chromatograph  analysis  time 

(ii)  The  process  time  delay 

(iii)  The  transport  time  delay  within  the 
sampling  system 

The  analysis  of  the  system  used  to  establish  the  total 
time  delay  is  given  in  Appendix  B. 


3 . 3  Compos i t i on  Measur emen t  s 

The  accuracy  of  terminal  composition  measurements  and 
the  reliability  of  the  measuring  devices  are  vital  to  the 
control  system. 

The  use  of  the  capacitance  probe  provides  continuous 
measurement  for  top  product  composition.  This  continuous 
signal  is  sampled  at  64  second  intervals  for  data  logging 
and  control  purposes.  This  method  of  analysis  provides  good 
accuracy  for  solutions  of  high  methanol  content  which  is  the 
case  for  the  top  product. 

A  Hewlett-Packard  model  5720A  gas  chromatograph  with  an 
automatic  sampling  valve  is  used  for  bottom  composition 
measurements  since  the  capacitance  method  of  analysis  is  not 
satisfactory  for  streams  with  a  low  methanol  content  (18). 

3.3.1  Process  and  Measurement  Noise 

High  frequency  noise  was  found  in  both  terminal 
compositions  as  can  be  seen  from  Figures  3.10  and  3.11.  The 
noise  is  in  part  due  to  the  measuring  devices,  but  it  is 
mainly  caused  by  the  noisy  manipulated  variables,  steam  and 
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Figure  3.10  Experimental  Steady  State  Top  Product 
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Figure  3.11  Experimental  Steady  State  Bot tom. Product 
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reflux  flows  as  shown  in  Figures  3.12  and  3.13. 

The  high  noise  level  in  the  steam  flow  rate  and  the 
bottom  composition  measurement  together  with  a  time  delay 
that  is  larger  than  the  system  time  constant  have  a 
detrimental  effect  on  the  performance  and  stability  of  the 
bottom  composition  control  system.. 

A  summary  of  all  the  measured  variables  of  the 
distillation  column  with  their  standard  deviation  is  given 
in  Table  3.4. 
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Figure  3.13  Experimental  Steady  State  Reflux  Flow  Noise 
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3.^  Summary  of  the  distillation  column  measured 

VARIABLES  (WITH  STANDARD  DEVIATION  BASED  ON 
50  SAMPLES  PER  VARIABLE  AT  STEADY  STATE 
CONDITIONS) 


FEED  FLOW 

18.012 

G/SEC 

DEV  =  0.0373 

REFLUX  FLOW  = 

10.304 

G/  SEC 

DE  V=  0.1301 

STEAM  FLOWW  = 

11.534 

G/SEC 

DEV  =  0.5275 

BOTTOM  PROD  = 

9.099 

G/  SEC 

DE  V=  0.2163 

TOP  PROD  = 

8.750 

G/SEC 

DEV=  0.0610 

COOL  WATER  =472.685 

G/SEC 

DE  V=  1.9260 

DIST  COMP 

96.583 

WT2 

ME  OH 

D  EV=  0.1  523 

BOTTOMS  C DM P= 

5.400 

WT£ 

ME  OH 

DE  V=  0.0000 

FEED  COMP 

49.250 

WT% 

MEOH 

DEV=  0.0000 

PRESSIR  E 

0.4  91 

K  PA 

DE  V  =  0.0585 

COND  LEVEL 

16.913 

CM 

DEV=  0.0449 

REB’R  LEVEL  = 

59.479 

CM 

DE  V=  0.3310 

DIFF  PRESS  = 

5.013 

KPA 

DEV  =  0.0041 

REB ’ R  0 ' HEAD= 

94.5 

DEG 

C 

DE  V=  0.214  5 

PLATE  1  TEMP  = 

85.3 

DEG 

c 

DEV=  0  . 2091 

PLATE  2  TEM  P= 

79.0 

DEG 

c 

DE V=  0.1661 

PLATE  3  TEMP  = 

75.  1 

DEG 

c 

DEV  =  0.1945 

PLATE  4  TEM P= 

74. 1 

DEG 

c 

DE  V  =  0.1881 

PLATE  5  T EMP  = 

70.  5 

DEG 

c 

DEV=  0.2120 

PLATE  6  TEM P= 

67.9 

DEG 

c 

DE  V=  0.2030 

PLATE  7  TEMP  = 

66.2 

DEG 

c 

DE V  =  0.1860 

PLATE  8  TEM P= 

64.5 

DEG 

c 

DE  V=  0.1861 

COMD  TEMP 

61.9 

DEG 

c 

DEV=  0.1794 

STEAM  TEMP  = 

89.7 

DEG 

c 

DE  V=  0.2476 

COMD  *  T  TEMP  = 

104.  4 

DEG 

c 

DEV=  0.2887 

REFLUX  TEMP  = 

48.5 

DEG 

c 

DE  V=  0.1984 

FEED  TEMP 

40.7 

DEG 

c 

DEV=  0.2053 

BOTTOMS  TEM P= 

50.9 

DEG 

c 

DE  V=  0.2142 

REB’R  TEMP  = 

93.  9 

DEG 

c 

DEV=  0.2300 

FEED  I  NLET  = 

72.9 

DEG 

c 

DEV=  0.2535 

REFLUX  INLET  = 

64.4 

DEG 

c 

OEV=  0.7267 

COL  O' HEAD  = 

64.9 

DEG 

c 

DEV=  0.1798 

WATER  INLET  = 

9.4 

DEG 

c 

DEV=  0.1782 

WATER  OUTL  ET  = 

21.3 

DEG 

c 

DE  V=  0.2021 

. 


4.  TOP  PRODUCT  COMPOSITION  CONTROL 


4.1  System  Description 

A  binary  solution  of  methanol  and  water  was  fed  to  the 
fourth  tray  of  the  eight  tray  distillation  column.  The  high 
methanol  concentration  overhead  vapor  was  condensed  in  the 
total  condenser.  Condenser  level  was  controlled  by  top 
product  flow  with  the  reflux  flow  chosen  to  control  the  top 
product  composition. 

Most  of  the  column's  instrumentation  is  pneumatic  with 
analog  controllers  in  each  of  the  flow  control  loops.  The 
column  was  operated  in  an  analog  supervisory  control  inode  as 
interfaced  to  the  IBM  1800  process  computer. 

4.2  Experimental  Evaluation 

4.2.1  System  Configuration  and  Implementation 

All  process  measurement  signals  were  transduced  to  the 
computer  for  logging  and/or  control  purposes.  The  control 
law  calculation  was  scheduled  to  run  every  64  seconds  with 
the  required  inputs  being  the  top  product  composition  and 
reflux  flow  rate.  The  IBM  1800  Direct  Digital  Control  (DDC) 
system  has  a  basic  scan  frequency  of  one  second.  Each  point 
in  the  data  base  was  scanned  at  the  specified  poll  time 
which  could  be  any  integer  multiple  of  the  basic  scan 
frequency.  The  raw  data  information  was  processed, 
converted,  and  put  in  the  data  base.  At  the  scheduled 
interval,  the  user  written  control  law  algorithm  calculated 
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the  error  difference  between  the  target  composition  set 
point  and  the  measured  composition  signal  and  then 
calculated  the  required  control  action  of  the  manipulated 
variable,  reflux  flow  rate.  Limit  checking  was  performed 
before  the  calculated  value  was  passed  for  transmission  as 
the  set  point  to  the  reflux  flow  controller.  Figure  4.1  is 
a  schematic  representat ion  of  the  system  which  shows  the 
interface  between  the  process  and  the  computer.  The  source 
listing  of  the  complete  set  of  control  programs  is  given  in 
Appendix  C. 

4.2.2  Control  Strategy 

4.2.2. 1  Proportional  Plus  Integral  Control 

In  order  to  achieve  a  valid  comparison  between  the  STR 
and  PI  control  performance,  the  control  performance  achieved 
using  a  well -tuned  PI  controller  was  used  as  a  reference  in 
evaluating  SIR  control  behavior.  The  controller’s  setting 
was  determined  using  the  Ziegler-Nichols  tuning  method  on 
the  actual  column.  Tuning  details  can  be  found  in  Appendix 
D. 

4. 2. 2. 2  Self-Tuning  Regulator 

With  the  real  time  data  base  and  the  information  vector 
of  previous  input  and  output  signals,  model  parameters  <^i 
and  ^i  were  estimated  at  every  sampling  interval.  At  the 
same  instant,  these  estimated  parameters  were  used  for 
calculation  of  the  control  signal  to  minimize  output 
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D/P  Differential  pressure  transmitter 

P/I  Pressure  to. current  converter 

I/P  Current  to  pressure  converter- 

PC  Analog  flow  controller 

SP  Set  point 

J'icAS  Measurement  " 

AT  Analyzer  Transmitter 

Figure  4.1  Process  &  Computer  Interface  for  Analog 

Supervisory  Control 
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variance  of  the  system.  Parameters  required  for  operation 
of  the  algorithm  that  were  selected  for  use  with  the  base 
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run  were: 


model  order 


N  =  2 


time  delay 

initial  parameter  values 
exponential  forgetting  factor 


scaling  factor 
initial  P ( 0)  value 


&o  =  0.01 

P ( 0 )  =  3,000  I 
K  =  1 
6(0)  =  0. 
f*.  =  1.0 


These  values  were  determined  by  digital  simulation 
prior  to  the  experimental  testing. 

Note  that  a  model  order  of  two  was  selected  to  allow 
for  variation  in  the  gain  and  time  constant  of  the  system 
despite  the  fact  that  a  first  order  plus  time  delay  model 
was  adequate  for  characterizing  a  particular  response.  This 
was  considered  necessary  because  the  gain  and  time  constants 
were  found  to  be  strongly  dependent  upon  the  direction  and 
magnitude  of  the  disturbances. 

4.2.3  Results 

Tests  were  performed  for  a  feed  flow  step  decrease  of 
twenty  percent  and  step  changes  in  composition  set  points  of 
+1.0%  and  -1.5%  from  the  steady  state  values.  The 
performance  of  both  the  PI  and  STR  control  algorithms  was 
evaluated  for  these  changes. 

When  the  column  reached  steady  state  at  the  required 
operating  condition,  the  control  task  was  initiated  via  the 
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teletype  terminal.  It  was  at  this  instant  that  run  time, 
time  and  magnitude  of  the  disturbance,  control  option,  and 
controller  constants  were  specified.  For  all  experimental 
runs,  a  run  time  of  180  time  units  was  used.  Each  time  unit 
spanned  64  seconds.  The  disturbance  or  change  in  set  point 
was  introduced  at  the  30th  time  unit  from  the  start  of  the 
run.  The  teletype  input  for  a  PI  control  test  run  for  a  1% 
set  point  increase  is  given  in  Figure  4.2. 

4.2.3. 1  Proportional  Plus  Integral  Control  Tests 

Ziegler -Nichols  tuning  of  the  top  product  composition 
control  yielded  a  control ler  setting  with  a  proportional 
gain  of  Kp=-45,  and  a  reset  time  T^=5.84  minutes.  A  summary 
of  the  excellent  results  is  given  in  Table  4.1  with  the 
controlled  responses  plotted  in  Figures  4.3  to  4.5. 

Table  4.1  Summary  Of  Top  Product  Experimental  PI 

Control  Tests 


Run 

Lr 

ISE 

F iqure 

Di s turbance 

E8 

-45. 

5.84 

0.0000 

4.3 

Feed  Flow  -20% 

E9 

-45. 

5.84 

0.0002 

4.4 

Set  Point  +  1% 

E 1 0 

-45. 

5.84 

0.0014 

4.5 

Set  Point  -  1 . 

4. 2. 3. 2  Experimental  Self-Tuning  Regulator  Tests 
i )  Effect  of  : 

Earlier  studi  es  (5)  had  indicated  that  the  choice  of 
was  not  crucial,  though  it  should  be  selected  in  a  manner  to 
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TIME  CMIN) 


TIME  (MIN) 


Figure  4.3  Experimental  PI  Control  Of  Top  Product 

Composition  For  a  -k0%  Step  Change  In  Feed  Flow 
Rate 


REFLUX  (G/S)  TCP  COMP  (PERCENT) 
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Figure  4.4  Experimental  P3  Control  Of  Top  Product 

Composition  ror  a  +1%  Step  Change  In  Set  Point 


REFLUX  (G/S)  TCP  CCMP  (PERCENT) 
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Figure  4.5  Experimental  PI  Control  Of  Top  Product 

Composition  For  a  -1.5%  Step  Change  In  Set  Point 
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avoid  instability.  However,  the  results  of  this  study 
indicate  that  the  choice  of  the  value  significantly 
influences  the  control  behavior.  The  tests  listed  in  Table 
4.2  were  made  to  study  the  effect  of  on  control 
performance.  The  output  responses  are  plotted  in  Figures 
4.6  to  4.11. 

Table  4.2  Effect  Of  fo  Value  On  Experimental  STR 

Control  Of  Top  Product  Composition 


M  =  2 


K  =  1 


L  =  3 


6(0)  =  0.0 


Run  |o 

Pj  0) 

ISE 

F igure 

Feed  Flow  -20% 

Ell  0.05 

-4 

250  I 

33 . x 1 0  ^ 

4.6 

E 1 2  0.01 

250  I 

22.x  10"T  • 

4.7 

E 1 3  0.005 

250  I 

14.x  10”* 

4.8 

Set  Point  +1% 

E  1 4  0.05 

3,000  I 

11.x  10'^ 

4.9 

E 1 5  0.01 

3,000  I 

7 . xl 0  ; 

4.  10 

E 1 6  0.005 

3,000  I 

6.x  10’* 

4.11 

The  results  in  all  runs  showed  that  the  control 
behavior  was  significantly  improved  by  lowering  the 
value.  The  maximum  deviation  of  the  composition  from  the 
set  point  value  for  the  feed  flow  rate  changes,  decreased  as 
the  value  of  was  decreased.  The  maximum  deviation  of  the 
composition  decreased  from  1.3  to  0.8  weight  percent  with  ^>o 
changed  from  a  value  of  0.05  to  0.005  as  shown  in  Figures 
4.6  and  4.8.  In  the  case  of  the  set  point  changes,  the 
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F i gure 


4.6  Experimental  STR  Control  Of  Top  Product 

Composition  For  a  -20%  Step  Change  In  Feed  Flow 
Rate,  ^0=  0.05,  P  ( 0 )  =  250^. 
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Figure  4.7  Experimental  STR  Control  Of  Top  Product 

Composition  For  a  -20%  Step  Change  In  Feed  Flow 
Rate,  &  =  0.01  ,  P(0)  =  250  1. 
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Figure  4.8  Experimental  5TR  Control  Of  Top  Product 

Composition  For  a  -20%  Step  Change  In. Feed  Flow 
Rate,  £0  =  0.005,  P ( 0 )  =  250 
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Figure  4.9  ExDerimental  STR  Control  Of  Top  Product 

Composition  For  a  +1%  Step  Chanae  In  Set  Point, 
=  0.05,  EJO)  =  3 0 0 0 
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Figure  4.10  Experimental  STR  Control  Of  Top  Product 

Composition  For  a  +1%  Step  Change  In  Set 
=  0.01,  P ( 0 )  =  3000  I. 
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Figure  4.11  Experimental  STR  Control  Of  Top  Product 

Composition  For  a  +1%  Step  Change  in  Set  Point, 
0.005,  E(0)  =  3000  1  . 
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overshoot  of  0.3%  observed  in  Figure  4.9  was  totally 
eliminated  by  lowering  j£<>  from  0.05  to  0.005  as  seen  in 
Figure  4.11.  For  both  the  servo  and  regulatory  control 
tests  performed,  the  best  control  behavior  was  obtained 
using  the  lowest  value.  The  control  performance  was  not 
as  sluggish  for  the  smaller  value  of  and  gave  the 
smallest  integral  square  error  (ISE)  value. 

In  examining  the  results  listed  in  Table  4.2,  it  can  be 
seen  that  an  increase  in  the  1/p0  value  reduces  the  ISE  for 
all  runs.  This  behavior  is  not  surprising  since  1/j30  acts 
as  a  gain  term  for  the  control  input  calculation  (cf 
Equation  2.11)  in  the  STR  algorithm.  An  increase  of  the 
1/p0  value  has  prompted  the  parameters  to  converge  faster. 
The  estimated  parameters  for  the  change  in  feed  flow  rate 
shown  in  Figure  4.12,  took  46  time  units  after  the 
disturbance  was  introduced  to  reach  their  steady  state 
values  with  =  0.05,  and  40  time  units  with  po  =  0.005  in 
Figure  4.14.  Figures  4.12  to  4.14  are  included  to 
demonstrate  typical  parameter  convergence. 

However,  the  reflux  flow  and  top  product  composition 
displayed  the  most  oscillatory  behavior  for  the  lowest 
value  used,  giving  a  "bang  bang"  like  performance.  This 
excessive  control  behavior  was  believed  to  be  caused  by  the 
high  1/p0  gain  effect  that  induced  the  high  sensitivity  of 
the  regulator  to  the  noise  in  the  system. 

In  examining  the  output  responses  plotted  in  Figures 
4.6  to  4.11,  it  can  be  seen  that  steady  state  error  (offset) 
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F igure 


4.12  Estimated  Parameters  For  Experimental  STR 

Control  Of  Top  Product  Composition  For  a  -20% 
Step  Change  In  Feed  Flow  Rate  With  60-  0.05, 
EJO)  =  250  U  v 
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Figure  4.13  Estimated  Parameters  For  Experimental  STR 

Control  Of  Top  Product  Composition  For  a  -20% 
Step  Change  In  Feed  Flow  Rate  With  0O  =  0.01, 
P(0)  =  250  1.  ' 
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Figure  4.14  Estimated  Parameters  For  Experimental  STR 

Control  Of  Top  Product  Composition  For  a  -20% 

Step  Change  In  Feed  Flow  Rate  With  =  0.005, 
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resulted,  with  the  offset  increasing  with  the  lower  values 
of  value  used.  Consider  the  STR  algorithm  as  given  by 
Equat ion  2.12 

u(t)  =  [  Om  (t)y(t)  +  ...+  y  ( t-n+ 1 )  ] 

-^»(t)u(t-1)  -  ...  -  ( t )  u  ( t  -  L  ) 

in  which  u(t)  and  y(t)  are  defined  as  the  normalized 
deviation  variables.  At  steady  state,  the  values  of  y(t) 
and  u(t)  will  be  the  same  as  their  past  values.  Hence,  by 
assuming  the  system  has  reached  steady  state  at  the  end  of 
the  run,  this  equation  can  be  rearranged  into  the  form 

y  =  u  *  |>o  (1.0  +2^1) 

Now  as  value  approaches  -1.0,  the  final  offset,  y  will 
approach  zero.  Table  4.3  is  a  summary  of  the  final 
estimated  o(i  and  parameters  for  the  control  tests  shown 
in  Figures  4.6  to  4.11. 
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Table  4.3  Effect  Of  Value  On  Estimated 

Parameter  Values 


Run  F iqure  Final _ Estimated  Parameter  Values 


Feed 

Flow  - 

20% 

c Cl 

Ell 

0.05 

4.12 

-0.3571 

-0.3490 

-0.7061 

-0.3383 

-0.3193 

-0.3013 

-0.9589 

E  12 

0.01 

4.  13 

-0.3287 

-0.3243 

-0.6530 

-0. 1870 
-0. 1832 
-0. 1796 

-0.5498 

E  1 3 

0.005 

4.  14 

-0.2487 

-0.2452 

-0.4939 

-0.0765 

-0.0750 

-0.0736 

-0.2251 

Set 

Point  +1% 

E  1 4 

0.05 

- 

-0.6256 

-0.4824 

-1  .1080 

-0.5140 
-0.3044 
-0. 1499 

-0.9683 

E  15 

0.01 

-0.4999 

-0.2924 

-0.7923 

-0.3155 

-0.2757 

-0.2744 

-0.8656 

E  16 

0.005 

-0.4699 

-0.2600 

-0.7299 

-0.2058 

-0.2095 

-0.2269 

-0.6422 

These 

resu 1 ts 

show  that 

as  £0  i  s 

decreased , 

the 

di f f erence 

between 

the  value  and 

-1.0  i ncreases .  The 

values  changed  from  -0.9589  to  -0.2251  in  Runs  Ell  to  E13 
for  a  20%  step  decrease  in  feed  flow  rate  and  from  -0.9683 
to  -0.6422  in  Runs  E14  to  E16  for  a  1%  increase  in  set 


point.  In  the  case  of  l,  the  absolute  value  decreased  as 
was  lowered,  thus  also  leading  to  an  increase  in  the 
offset.  The  combined  effect  of  £cCc  and  outweighed  the 

effect  of  l^o  as  shown  by  the  larger  offset  that  resulted. 


- 

. 
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For  example,  in  the  case  of  Run  E12,  the  values  were  = 

-0.5498;  =  0.6530;  =  0.01;  giving 

y  =  (12  -  10.4)/(10.4  WO.On*  (1  -  0 . 5498 ) /-0 . 6530  =  0.00106 

so  for  the  set  point  of  96.2%,  the  final  value  of  y  wi 1 1  be 
96.1%.  This  calculated  value  is  in  agreement  with  the  final 
value  in  Figure  4.7.  This  trend  is  consistent  for  all  of 
the  tests. 

For  performance  evaluation,  since  an  ISE  index  was 
used,  it  should  be  noted  that  the  duration  of  the  run 
influences  the  accumulated  error  due  to  steady  state  offset. 
For  the  tests  performed,  one  could  achieve  faster  response 
and  reduced  ISE  by  lowering  at  the  expense  of  increasing 
offset . 

* 

ii)  Effect  of  £(0)  value: 

The  P(0)  value  set  at  the  beginning  of  the  run  was 
determined  by  the  amount  of  confidence  which  could  be  placed 
in  the  initial  parameter  estimates.  If  the  parameter 
estimates  were  considered  to  be  close  to  the  true  values, 
the  P(0)  value  was  set  to  a  low  value.  Tests  performed 
using  different  values  of  JP  ( 0 )  are  summarized  in  Table  4.4. 
The  output  responses  are  plotted  in  Figures  4.15  to  4.20. 
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Effect  Of  Initial  Value  Of  Covariance  Matrix, 
£(0)  On  The  Experimental  STR  Control  Of  Top 
Product  Composition 


Run 

|o 

JL(0) 

ISE 

Figure 

Set 

Point  +1% 

E  1 5 

0.01 

3,000  I 

7.3x1 oT 
6.0x107 

4.  10 

E  1 7 

0.01 

10,000  JL 

4.  15 

E  1 8 

0.01 

30,000  JL 

5. 1x1  CP 

4.  16 

Set 

Point  -1.5% 

E  1 9 
E20 

0.01 

0.01 

10,000  ± 
30,000  ± 

11.9x1  (P 
14.4x10^ 

4.  17 
4.18 

Feed 

Flow  -20% 

E  1 3 
E21 

0.005 

0.005 

250  I 
500  1 

1 4  .  Ox  1  0 
8.7x10j 
4.7x1 (j 

4.8 

4.19 

E22 

0.005 

1,000  I 

4.20 

< 
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F igure  4.15 


Experimental  STR  Control  Of  Top  Product 
Composition  For  a  +1%  Step  Change  In  Set 
=  0.01 ,  P(0)  =  10000  I. 
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Figure  4.16  Experimental  STR  Control  Of  Top  Product 

Composition  For  a  +1%  Step  Change  In  Set  Point, 
0o  =  0.01 ,  P(0)  =  30000  JL 
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Figure  4.18  Experimental  STR  Control  Of  Top  Product 

Composition  For  a  -1.5%  Step  Change  In  Set 
Point,  R0  =  0.01  ,  PJO)  =  30000  r. 


' 


REFLUX  (G/S)  TCP  COMP  (PERCENT) 


63 


TIME  (MIN) 


.Figure  4.19  Experimental  STR  Control  Of  Top  Product 

Composition  For  a  -20%  Step  Change  In  Feed  Flow 
Rate,  &>  =  0.005,  P(0)  =  500  ^ 
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Figure  4.20  Experimental  STR  Control  Of  Top  Product 

Composition  For  a  -20%  Step  Change  In  Feed  Flow 
Rate,  =  0.005,  PjO)  =  1000  K 
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From  the  results  obtained,  the  sustained  oscillations 
that  have  been  evident  in  all  STR  runs  became  more 
pronounced  when  a  higher  PJ  0 )  value  was  used.  This  is 
obvious  from  the  responses  plotted  in  Figures  4.15  -  4.20. 
When  P(0)  =  30,000  I_  was  used,  as  can  be  seen  in  Figure 
4.18,  the  manipulated  variable  (reflux)  banged  between  its 
high  and  low  limits  at  the  time  the  set  point  change  was 
introduced.  This  behavior  resulted  because  of  large 
fluctuations  in  the  estimated  parameter  values  as  shown  in 
Figure  4.21.  For  a  1%  step  increase  in  set  point,  with  PJ0) 
=  30,000  I_,  the  control  behavior  became  oscillatory  as  shown 
by  the  response  in  Figure  4.16.  The  poor  control  is  to  be 
expected  due  to  the  changes  in  the  parameter  values  (cf 
Figure  4.22).  The  parameters  drift  away  from  each  other, 
unlike  the  adaption  pattern  of  the  parameters  of  the  other 
tests,  and  never  reach  steady  values. 

Changing  the  value  of  P(0)  showed  that  the  effect  on 
controller  performance,  to  a  large  extent,  was  almost 
identical  to  the  gain  term  in  a  proportional  controller.  In 
addition,  the  steady  state  offset  was  found  to  decrease  as 
the  £(0)  value  was  increased.  From  the  and  values 
given  in  Table  4.5,  it  can  be  seen  that  as  the  value  of  £(0) 
increased,  £-^2  moved  closer  to  -1.0  and  the  absolute  value 
of  i-oLl  increased,  so  for  the  same  value,  the  offset  would 
certainly  be  smaller.  The  1SE  values  were  lowered  by 
increasing  JP  ( 0 )  for  all  tests  except  the  1.5%  step  decrease 
in  set  point  run  which  was  caused  mainly  by  the  initial 
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Figure  4.21  Estimated  Parameters  For  Experimental  STR 

Control  Of  Top  Product  Composition  For  a  -1.5% 
Step  Change  In  Set  Point,  =  0.01,  JPJ  0 )  = 
3000  I  .  ' 
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Figure  4.22  Estimated  Parameters  For  Experimental  STR 

Control  Of  Top  Product  Composition  For  a  +1% 
Step  Change  In  Set  Point,  p0=  0.01,  P_(  0 )  = 
30000  I. 
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Table  4.5  Effect  Of  E.(0)  Value  On  The  Estimated 

Parameters 


Run 

PjO) 

F i gure 

Final 

Estimated 

Parameter 

V  a  1  ue 

q<i 

& 

Set 

Point  +1% 

\ 

E  1 5 

3,000  I 

4.10 

-0.4999 

-0.2924 

-0.7923 

-0.3155 

-0.2757 

-0.2744 

-0.8656 

E  1 7 

. 10,000  I 

-0.6945 
-0. 1668 

-0.8613 

-0.3454 

-0.2969 

-0.2977 

-0.9400 

E  1 8 

30,000  l 

4.22 

-0.9795 

0.1161 

-0.8634 

-0.5376 

-0.1979 

-0.2100 

-0.9464 

Set 

Poi nt  -1.5% 

E 1 9 

10,000  i 

- 

-0.8636 

-0.0312 

-0.8324 

-0.2964 

-0.2439 

-0.2973 

-0.8376 

E20 

30,000 

4.21 

-1.0953 

0.2205 

-0.8748 

-0.0348 

-0.3866 

-0.4283 

-0.8497 

Feed 

Flow  -20% 

• 

E 1 3 

250  l 

4.8 

-0.2487 

-0.2452 

-0.4939 

-0.0765 

-0.0750 

-0.0736 

-0.2251 

E21 

500  I 

a 

— 

-0.2899 

-0.2834 

-0.5733 

-0.  1132 
-0. 1083 
-0.1083 

-0.3298 

E22 

1,000  1 

-  . 

-0.2980 

-0.2888 

-0.5868 

-0.1213 
-0. 1185 
-0. 1170 

-0.3568 

' 
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4.3  Simulation  Study 

4.3.1  System  Configuration  and  Implementation 

The  simulations  were  performed  using  the  University  of 
Alberta  computing  services  Amdahl  470V/6  computer  system. 

The  nonlinear  distillation  process  model  developed  by 
Simonsmeier  (20)  was  incorporated  into  the  simulation 
program  to  generate  the  distillation  column  dynamic 
behavior.  The  PI  and  STR  control  routines  are  identical  to 
those  used  in  the  experimental  study  with  implementation  of 
the  control  action  on  the  model  performed  in  a  manner 
analogous  to  that  of  the  experimental  study.  The  source 
listing  of  the  complete  set  of  control  programs  used  in  the 
simulations  is  given  in  Appendix  E. 

Prior  to  any  simulation  run,  an  input  data  file  was  set 
up  to  provide  the  initial  steady  state  operating  conditions 
for  the  column  and  the  control  data  parameters  required.  A 
run  time  of  169  time  units,  with  each  time  unit  being  64 
seconds,  was  used  for  all  the  simulations.  The  disturbance 
was  always  introduced  after  20  time  units  of  simulation  time 
had  elapsed.  The  program  integration  interval  was  fixed  at 
8  seconds.  The  integration  interval  was  selected  as  a 
compromise  between  reducing  the  computer  execution  time  and 
accurately  calculating  the  transient  response  of  the  column. 
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4.3.2  Proportional  Plus  Integral  Control  Tests 

The  Ziegler -Nichol s  tuning  constants  of  Kp  =  -45.,  T^  = 
5.84  min.  determined  exper imenta 1 1 y ,  were  used  for  the 
simulation  runs.  The  results  are  summarized  in  Table  4.6. 

Table  4.6  Summary  Of  Simulation  Results  For  Top  Product 

Composition  PI  Control 


Run 

Di sturbance 

i<P 

Ir 

ISE 

F iqure 

S3 

Feed  Flow  -20% 

-45. 

5.84 

1 

.0x1 0"^" 

4.23 

S4 

Set  Point  +1.0% 

-45. 

5.84 

21 

.oxio'j 

4.24 

S5 

Set  Point  -1.5% 

-45. 

5.84 

59 

.0x1  O'4 

4.24 

Operating  Conditions 

Feed  Flow  18.00  g/s  Top  Composition  95.7% 

Reflux  Flow  11.00  g/s  Bottom  Composition  5.6% 

Steam  Flow  11.95  g/s 

For  a  20%  step  decrease  in  feed  flow  rate,  the  PI  controller 
gave  excellent  control  behavior  as  can  be  seen  from  the 
response  in  Figure  4.23.  The  maximum  deviation  was  only 
0.07  wt  %  from  the  set  point  value  so  further  improvement  in 
the  control  performance  by  further  tuning  will  likely  not 
result.  For  step  changes  in  set  point,  satisfactory  control 
responses  were  also  obtained.  As  can  be  seen  from  the 
control  behavior  in  Figures  4.24  and  4.25,  that  resulted  by 
using  the  same  controller  settings  as  for  the  disturbance  in 
feed  flow  rate,  a  larger  overshoot  was  observed  for  the 
decrease  in  set  point  as  compared  to  the  increase.  This  is 
due  to  the  nonlinear  dynamic  character i st ics  of  the  column. 
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Figure  4.23  Simulation  Of  PI  Top  Product  Composition 

For  a  -20%  Step  Change  In  Feed  Flow  Rate 
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Figure  4.24 


Simulation  Of  PI  Top  Product  Composition 
For  a  +1%  Step  Change  In  Set  Point 
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Figure  4.25  Simulation  Of  PI  Top  Product  Composition  Control 

For  a  -1.5%  Step  Change  In  Set  Point 
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4.3.3  Self-Tuning  Regulator  Tests 

Preliminary  simulations  were  carried  out  to  determine  a 
suitable  set  of  numerical  values  for  the  STR  parameters. 

The  operating  conditions  were  set  identical  to  those  used  in 
an  earlier  study  by  Sastry  et  al  (15).  The  following 
parameters  values  as  used  by  Sastry,  were  found  to  be 
satisfactory. 

=  0.01 
jU.  =  1.00 
P( 0)  =  3,000  l 
0(0)  =  0.0 

The  output  responses  for  the  STR  runs  using  this  set  of 
parameters  are  plotted  in  Figures  4.26  to  4.28,  and  the 
results  tabulated  in  Table  4.7. 
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Table  4.7  Summary  Of  Simulated  STR  Top  Product 

Composition  Control  Results 


M 

=  2 

k 

=  1 

L  =  2 

=  0.01 

P(0)  =  3 

,000  I. 

Run  ISE 

Final 

XD 

Final  Par ameter s 

F i gure 

Di sturbance 

-4 

OCL 

S8  108.6x10 

95. 

95% 

-0.4333 

-0.2410 

4.26 

Feed  Flow 

-0.4055 

-0.2318 

-20% 

1-0.8388 

1 0.4728 

S9  33.7x10 

96. 

95% 

-0.4529 

-0.0739 

4.27 

Set  Point 

-0.2356 

-0.0592 

+  1% 

1-0.6885 

1-0.  1331 

-4 

S11  69.1x10 

94. 

57% 

-0.5680 

-0.0816 

4.28 

Set  Point 

-0.  1060 

-0. 1216 

-1.5% 

1-0.6740 

1-0.2032 

Operatinq  Conditions 

Feed  Flow 

18.00 

g/s 

Top 

Compos i t i on 

=  96.02% 

Reflux  Flow 

17.90 

g/s 

Bottom  Compositiion 

=  1.49% 

Steam  Flow 

15.40 

g/s 

During  the  initial 

period  of 

operation,  the 

parameters 

for  a  1.5%  step  decrease  in  set  point  control  exhibited  a 


large  fluctuation  causing  the  manipulated  variable  (reflux 
flow)  to  be  driven  to  its  control  limit  at  6.0  g/s  before 
levelling  off  to  a  new  steady  state  level.  This  behavior 
was  similar  to  that  observed  in  the  subsequent  experimental 
work.  The  control  performance  was  reasonably  satisfactory, 
considering  zero  initial  parameter  estimates  were  used, 
except  that  offsets  were  observed  in  the  product  composition 
for  either  step  changes  in  load  or  set  point, as  was  the  case 
with  the  experimental  results. 


Reflux  (g/s)  Top  Comp  (.%) 
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Figure  4.26  Simulation  Of  STR  Top  Product  Composition 

Control  For  a  -20%  Step  Change  In. Feed  Flow  Rate 
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Figure  4.28  Simulation  Of  STR  Top  Product  Composition 

Control  For  a  - 1 . 5%  Step  Change  In  Set  Point 
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In  the  search  for  a  P ( 0 )  value  that  could  be  used  for 
both  load  disturbances  and  set  point  changes,  tests  were 
performed  to  study  its  effect  and  the  results  are  summarized 
in  Table  4.8.  It  was  found  that  the  P(0)  value  could  be 
increased  to  40,000  I  for  the  set  point  changes  without 
causing  unstable  responses.  However  for  the  feed  flow  rate 
decrease  of  -20%,  a  value  of  P(0)=3,000  £  appeared  to  be  the 
upper  limit  to  avoid  oscillatory  behavior. 
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The  simulation  results  indicated  that  the  effects  of 
the  P.(0)  value  on  offset  and  integral  square  error  were 
consistent,  with  that  found  in  the  experimental  runs.  Both 
the  ISE  value  and  offset  decreased  with  a  larger  JP  ( 0 )  value. 
The  effect  on  offset  could  again  be  explained  by  the  fact 
that  for  an  increase  of  P.(0)  from  3,000  (Run  S9)  to  40,000  J. 
(Run  S10)  that  had  increased  from  -0.1231  to  -0.5306 

and  from  -0.6885  to  -0.9069. 

Another  STR  parameter,  the  number  of  ^  parameters, 

(i.e.  the  L  value)  was  evaluated.  Use  of  a  value  of  3 
instead  of  2  was  investigated,  but  only  minor  improvement 
resulted  as  can  be  seen  from  the  results  in  Table  4.9. 
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Table  4.9 


Summary  Of  The  Effect  Of  The  Number  Of  ?> 
Parameters  (L)  On  The  Simulated  STR  Control  „ 
Of  Top  Product  Composition 


IV!  =  2 


K  =  1 


o 

II 

0.01 

PJ0)  = 

40, OOOl 

Run 

L 

Final  XD 

ISE 

Final  Parameters 

Set 

Point  +1 .0% 

o(; 

*oU 

Jk 

\ 

0.3089 

S10 

2 

96.98% 

4 

26.0x10 

-0.8745 

-0.9069  - 

-0.5306 

-0.0324 

- 

0.2217 

S 1  3 

2 

96.99% 

A 

25.6x10 

-0.8745 

-0.8608  - 

0.2541 

-0.7203 

-0.0098 

- 

0.2188 

- 

0.2374 

Set 

Point  -1.5% 

S12 

2 

94.52% 

-4 

62.0x10 

-0.9994 

-0.7557  - 

0.3342 

-0.7440 

+0.2437 

- 

0.4098 

S 1 4 

3 

94.51% 

61.3x10* 

-0.9686 

-0.6987  - 

0.2415 

-0.8861 

+0.2699 

- 

0.3284 

-0.2662 
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The  ISE  and  offset  were  reduced  somewhat,  but  the  final 


estimated  values  of 


parameters  adapted  in  such 


a  fashion  that  the  effect  on  offset  of  their  summation 
values  cancel  each  other.  These  results  suggest  that  when  a 
higher  than  adequate  number  of  ^  parameters  is  used,  only 
marginal  improvement  may  result  and  then  at  the  expense  of 
additional  computation  effort. 

4.4  Discussion 

The  top  product  control  system  is  fast  and  stable.  It 
has  a  short  time  delay  and  is  relatively  noise  free.  The 
operating  environment  is  unchanging.  These  combined  factors 
allow  for  good  control  behavior  using  a  conventional 
controller  so  i t  is  not  surprising  that  there  is  no 
advantage  in  using  an  advanced  control  scheme  like  the 
self-tuning  regulator  which  is  mainly  designed  to  deal  with 
processes  that  have  a  time  delay,  are  subject  to  a  changing 
environment,  stochastic  disturbances  and  system  noise. 

For  control  of  top  product  composition,  the 
conventional  proportional  plus  integral  controller  once  well 
tuned  performed  superbly  in  both  simulation  and  experimental 
tests.  The  responses  for  a  -1.5%  step  change  in  set  point 
as  shown  in  Figures  4.5  and  4.25  may  have  to  be  detuned  to 
reduce  the  overshoot,  but  overall  PI  control  provides  a 
control  performance  that  is  satisfactory.  In  both 
simulation  and  experimental  tests,  PI  clearly  outperformed 
STR,  yet  in  a  previous  study  by  Sastry  el  al  (15),  in  a 
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comparison  between  the  SIR  and  a  moderately  tuned  PI 
controller,  it  was  found  that  the  SIR  provided  the  better 
performance.  It  should  be  noted  that  a  comparison  between 
well  tuned  PI  control  and  STR  control,  with  zero  initial 
parameter  estimates,  is  somewhat  unfair  since  in  general  the 
STR  begins  with  very  poor  initial  parameter  estimates  while 
the  PI  controller  has  previously  been  tuned,  usually  for  the 
disturbance  of  interest.  In  addition,  using  the  integral 
square  error  as  the  performance  index  gives  a  severe  measure 
of  STR  performance  since  the  initial  parameter  estimation 
period  gives  rise  to  a  relatively  high  integral  square  error 
before  the  model  acquires  a  reasonable  set  of  parameters. 

The  results  have  shown  that  the  choice  of  and  the 
P(0)  value  significantly  influence  control  performance. 
Pretuning  of  these  parameters  would  be  advantageous  before 
using  the  STR  for  control  testing.  The  bang  bang  type  of 
control  performance  that  prevails  with  a  high  f^(0)  value 
and/or  a  small  value  suggests  that  by  adding  control 
weighting  as  proposed  by  Gawthrop  and  Clarke  (17)  the 
control  performance  should  improve.  The  problem  of  offset 
could  be  resolved  by  increasing  the  £(0)  value,  or  replacing 
u(t)  by  u(t)  -  u ( t - 1 )  in  equations  (2.10  and  2.11)  (5). 
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5.  BOTTOM  PRODUCT  COMPOSITION  CONTROL 


5.1  Description  of  Simulation  Study 

The  column  model  used  in  the  simulation  program  was 
identical  to  the  one  used  in  the  top  product  control  study 
except  for  a  program  counter  incorporated  into  the  program. 
This  counter  delayed  the  bottom  composition  calculated  by 
the  model  for  4  time  units  before  the  control  routine  was 
allowed  access  to  the  information.  For  the  bottom  product 
composition  control,  a  simulation  run  time  of  45  time  units 
was  employed  with  each  time  unit  being  256  seconds. 

The  counter  was  implemented  in  the  program  in  order  to 
simulate  the  equivalent  time  delay  present  in  the  bottom 
product  G.C.  analysis  system.  The  effect  of  changes  in 
steam  flow  were  not  delayed  but  were  input  directly  to  the 
model  to  immediately  influence  the  bottom  composition.  The 
control  of  the  bottom  dynamic  behavior  of  the  column  was 
complicated  by  this  excessive  transport  time  delay. 

The  control  law  calculation  was  scheduled  to  run  every 
256  seconds  with  the  required  inputs  being  the  bottom 
product  composition  and  steam  flow  rate.  The  operating 
sequence  of  events  in  the  cycle  can  be  summarized  as 
follows: 

1.  Bottom  composition  is  updated  by  the  simulated  gas 
chromatograph  analysis. 

2.  Values  for  steam  flow  and  bottom  conposition  are 
acquired  by  the  control  program. 
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3.  The  control  signal  is  calculated  using  the  specified 
control  routine. 

4.  The  control  signal  is  cascaded  to  the  steam  flow  analog 
control ler . 

5.2  Proportional  Plus  Integral  Control 

The  initial  controller  settings  were  determined 
experimentally  using  the  Ziegler-Nichols  tuning  method. 
Details  of  the  tests  performed  are  given  in  Appendix  F.  The 
constants  were  further  tuned  for  a  20%  step  decrease  in  feed 
flow  rate.  A  summary  of  the  different  PI  control  tests  is 
provided  in  Table  5.1. 
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Table  5.1  Summary  Of  Simulation  Results  For  PI 

Control  Of  Bottom  Product  Composition 


Run 

!<P 

V 

ISE 

B 1 

0.00 

0. 

3.297 

B2 

0.06 

42. 

1.596 

B3 

0.08 

42. 

1 .  167 

B4 

0.06 

42. 

1.925 

B5 

0.08 

42. 

1.575 

B6 

0.06 

42. 

0.082 

B7 

0.08 

42. 

0.081 

B8 

0.06 

42. 

0.082 

B9 

0.08 

42. 

0.090 

Operat 

inq  Conditions 

Feed  F 

low 

18.00 

g/s 

Reflux 

F  low 

11.00 

g/s 

Steam 

F  low 

11.95 

g/s 

Figure  Pi sturbance 


5.1 

Feed  Flow 

-20% 

5.1 

Feed  Flow 

-20% 

5.  1 

Feed  Flow 

-20% 

5.2 

Feed  Flow 

+20% 

5.2 

Feed  Flow 

+20% 

5.3 

Set  Point 

-  2% 

5.3 

Set  Point 

-  2% 

5.4 

Set  Point 

5^ 

CM 

+ 

5.4 

Set  Point 

+  2% 

Top  Composition  =  95.7% 
Bottom  Composition  =  5.6% 
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5.2.1  Load  Disturbance  Control 

All  feed  flow  rate  disturbances  were  introduced  2.5 
time  units  from  the  start  of  the  simulation  run.  For  a  20% 
step  decrease  in  feed  flow  rate  as  shown  in  Figure  5.1, the 
bottom  composition  did  not  start  decreasing  until  4  time 
units  after  the  disturbance  was  introduced  due  to  the  time 
delay  incorporated  into  the  output  from  the  column  dynamic 
model.  It  was  essentially  after  this  elapsed  time  delay 
that  the  manipulated  variable  (steam  flow)  was  gradually 
lowered  to  raise  the  composition  back  to  its  set  point.  The 
output  responses  were  sluggish,  due  to  the  long  system  time 
delay  and  time  constant.  The  open  loop  response  has  been 
included  to  provide  a  contrast  to  the  control  runs. 

Clearly,  of  the  different  simulations  per formed , the 
preferred  controller  setting  was  Kp  =  0.08,  T^  =  42  minutes 
as  shown  by  the  lower  ISE  value  of  1.167  as  listed  in  Table 
5. 1 . 

For  a  feed  flow  rate  step  increase  of  20%,  this  same 
controller  setting  again  provided  the  better  control 
performance  as  shown  by  the  responses  plotted  in  Figure  5.2. 
As  can  be  seen,  the  system  is  easier  to  control  for  an 
increase  in  feed  flow  rate  than  a  decrease  in  feed  flow 
rate.  This  is  due  to  the  different  gains  and  time  constants 
that  exist  due  to  the  nonlinear  dynamic  response  of  the 
column  to  positive  and  negative  feed  flow  rate  disturbances. 
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Figure  5.2  Simulated  Behavior  For  PI  Control  Of  Bottom 

Product  Composition  For  a  20%  Step  Increase  In 
Feed  Flow  Rate 


■  ■  . 

. 


92 


5.2.2  Set  Point  Control 

Servo  control  performance  was  evaluated  for  a  -2%  step 
change  in  set  point  with  the  resulting  manipulated  and 
controlled  variable  responses,  for  the  two  selected 
settings,  given  in  Figure  5.3.  As  can  be  seen,  both  choices 
of  settings  provided  reasonably  satisfactory  control  with 
almost  identical  ISE  values  as  shown  in  Table  5.1. 

For  a  +2%  step  change  in  set  point,  the  control 
performance  using  Kp  =  0.08,  =42.  resulted  in  an 

overshoot  and  a  slightly  oscillating  behavior.  By  lowering 
Kp  to  the  value  of  0.06,  an  improved  control  performance  was 
obtained  (cf  Figure  5.4). 

On  the  basis  of  the  simulation  results  listed  in  Table 
5.1,  the  control  performance  obtained  using  the  controller 
setting  of  Kp  =  0.08,  TR  =  42.  will  be  used  as  the  base  case 
for  evaluation  of  the  self-tuning  regulator  for  both  servo 
and  regulatory  control. 


5.3  Self-Tuning  Regulator 

5.3.1  System  Identification 

Preliminary  simulations  were  carried  out  for  control  of 
bottom  composition  using  the  STR  with  zero  initial 
estimates,  but  difficulties  were  encountered  in  trying  to 
obtain  a  converged  set  of  parameter  estimates.  This  is 
illustrated  by  the  results  of  Run  B10.  For  a  feed  flow 
disturbance  of  -20%,  using  =  10.0,  P_(  0 )  =  50.1  and  0(0/  = 
0.0  for  a  total  run  time  of  700  minutes,  the  responses  and 
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parameter  behavior  are  given  in  Figures  5.5  and  5.6 
respectively.  As  can  be  seen, the  system  tends  to  go 
unstable  with  the  bottom  composition  response  oscillating 
with  an  increased  magnitude. 

For  the  model  chosen  with  seven  parameters,  five  and 

two  oCi  ,  must  be  estimated.  As  Figure  5.6  shows,  the 
estimated  parameters  were  diverging  from  one  another, 
cycling  at  the  same  frequency  as  the  oscillations  found  with 
the  steam  flow  and  bottom  composition  in  Figure  5.5.  Since 
the  parameters  exhibited  large  initial  fluctuations  and 
then  levelled  off  to  some  steady  state  values  gradually,  the 
transient  behavior  evident  in  Figure  5.5  must  be  attributed 
to  the  behavior  of  o/;  estimates. 

This  problem  in  system  identification  i-s  believed  to  be 
caused  by  the  large  time  delay  relative  to  the  time 
constant.  Instead  of  using  zero  initial  parameter  estimates 
that  failed  to  yield  good  control,  a  multiple  pass  method 
was  used.  By  performing  parameter  estimation  runs  prior  to 
the  actual  control  runs,  it  is  expected  that  better 
performance  due  to  a  better  set  of  initial  estimated 
parameters,  will  result.  For  Run  B 1 1 ,  a  square  pulse 
disturbance  in  feed  flow  rate  was  introduced  into  the 
system.  This  is  accomplished  by  introducing  a  20%  step 
decrease  in  feed  flow  rate  after  10  minutes  of  simulation 
time  had  elapsed.  The  decrease  in  feed  flow  rate  was  then 
maintained  for  a  duration  of  170  minutes.  The  system 
response  and  the  parameters  are  plotted  in  Figures  5./  and 
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Figure  5.5  Simulated  Response  For  STR  Bottom  Product 

Composition  Control  For  a  -20%  Step  Change  In 
Feed  Flow  Rate , Star t ing  With  Zero  Initial 
Parameter  F.stimates 
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Figure  5.6  Parameter  Behavior  For  Simulated  STR  Bottom 

Product  Composition  Control  For  a'  -20%  Step 
Change  In  Feed  Flow  Rate,  Starting  With  Zero 
Initial  Parameter  Estimates 
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5.8  respectively.  The  2nd  arrow  indicates  the  time  that  the 
flow  is  returned  to  its  initial  steady  state  level.  As 
would  be  expected,  the  system  response  was  identical  to  that 
of  Run  BIO  until  the  feed  was  returned  to  its  normal  level 
after  180  minutes  of  simulation  elapsed  time.  The  behavior 
of  the  parameters  (cf  Figure  5.8)  shows  that  as  the  feed  is 
returned  to  its  normal  value,  the  04  parameters  start 
converging  instead  of  diverging  as  in  Run  BIO.  Despite  this 
convergence,  the  bottom  product  composition  response 
exhibited  significant  oscillation  but  did  manage  to  return 
to  its  set  point  after  about  700  minutes  of  time. 

By  performing  the  parameter  estimation  run  (Run  B 1 1 ) ,  a 
set  of  converged  estimated  parameters  is  obtained. 
Furthermore,  the  fact  that  the  output  returned  to  its  set 
point  provides  some  confidence  in  the  estimated  parameter 
values.  Results  from  parameter  estimation  runs  performed 
using  different  flow  patterns  are  discussed  in  Appendix  G. 

With  the  estimated  parameters  obtained  from  using  the 
square  pulse  disturbance  in  Run  B 1 1 ,  the  parameters  were 
further  adapted  using  a  20%  step  decrease  in  feed  flow  rate 
in  Run  B12.  The  results  in  Figure  5.9  show  that  the  bottom 
composition  started  to  move  back  to  its  set  point  as  the 
estimated  parameters  began  to  reach  their  steady  state 
values  in  180  minutes  or  45  control  intervals  after  the 
disturbance  was  introduced.  This  is  a  significant 
improvement  compared  to  Run  B10  in  which  use  of  the  zero 
initial  parameter  estimates  resulted  in  an  unstable  system. 
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Figure  5.7  Simulated  Response  For  STR  Bottom  Product 

Composition  Control  For  a  Pulse  Disturbance  In 
Feed  Flow  Rate  Of  '3.6  g/s  In  Magnitude  For  a 
Duration  Of  170  Minutes,  Starting  With  Zero 
Initial  Parameter  Estimates 
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Figure  5.8  Adaption  Pattern  Of  Parameters  For  a  Pulse 

Disturbance  In  Feed  Flow  Rate  Of  -3.6  g/s  In 
Magnitude  For  a  Duration  Of  170  Minutes,  Starting 
With  Zero  Initial  Parameter  Estimates 
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The  parameters  identified  in  Run  B12  were  such  that  the 
^p'Lwas  close  to  the  desired  value  of  -1.0  which  implies 
that  reliable  parameter  est imates. have  been  obtained.  The 
behavior  of  the  parameter  estimates  in  Run  B12  are  plotted 
in  Figure  5.10  and  the  results  for  all  the  parameter 
estimation  runs  are  summarized  in  Table  5.2. 

Table  5.2  Summary  Of  STR  Parameter  Identification 

Runs  For  Bottom  Product  Composition  Control 


M  =  2  K  =  4  L  =  5 


_E(0)  =  50. 


-Bun 

F i aure 

Initial 

final 

F  low 

Par ameter s 

Parameters 

Di sturbance 

Sl 

B10 

5.5 

0.0 

-5.0358 

-0.4895 

18.00 

g/s 

5.6 

3.3362 

0.1190 

0.2624 

o- 

10 

700  min. 

-0. 1594 

14.40 

g/s 

0.0596 

B 1 1 

5.7 

0.0 

-0.6066 

0.0398 

18.00 

g/s 

5.8 

-0.8923 

0.3613 

0.3830 

o- 

10 

180 

700  min. 

' 

0.0480 

0.6255 

14.40 

g/s 

B 12 

5.9 

Final 

-4.1887 

0. 1503 

18.00  g/s 

5.10 

Parameters 

3.0602 

0.0120 

or 

Of  Run  B 1 1 

0.1294 
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700  min. 

-0.3217 
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-0.959515 
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Figure  5.9  Simulated  Response  Of  STR  Bottom  Product 

Composition  Control  For  a  -20%  Step  Change  In 
Feed  Flow  Rate,  Starting  With  The  Initial 
Estimates  As  Those  Obtained  From  Run  B11 
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Figure  5.10  Parameter  Behavior  For  Simulated  STR  Bottom 

Composition  Control  For  a  -20%  Step  Change  In 
Feed  Flow  Rate,  Starting  With  The  Initial 
Estimates  As  Those  Obtained  From  Run  B11 
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5.3.2  Load  Disturbance  Control 

With  the  final  estimated  parameters  obtained  from  Run 
B 1 2 ,  the  control  performance  for  a  load  disturbance  was 
tested  by  introducing  a  20%  step  decrease  in  feed  flow  rate. 
In  these  tests,  the  parameter  identi f ication  was  stopped 
prior  to  the  introduction  of  the  disturbance.  This 
approach .designated  as  Fixed  Parameter  Control  (FPC),  has 
resulted  in  good  control  performance  with  the  output 
response  significantly  better  than  the  well-tuned  PI.  To 
compare  the  performance  of  the  two  types  of  control ,  the  PI 
and  FPC  responses  are  plotted  in  Figure  5.11.  The  FPC 
response  is  quicker  and  the  resulting  ISE  is  lower  compared 
to  the  value  obtained  under  PI  control.  The  results  for  the 
step  changes  in  feed  flow  rate  tests  are  summarized  in  Table 


5.3. 

Table  5.3  Summary  Of  Multiple  Pass  FPC  Control 


Of 

In 

Bottom  Composition  For  Step  Changes 

Feed  Flow  Rate 

Run 

ISE 

Figure 

Di s turbance 

B13 

0.792 

5.11 

Feed  Flow  -20% 

B 14 

0.835 

5.12 

Feed  Flow  +20% 

For  a  +20% 

feed  flow 

rate  load  disturbance,  with  the 

same  set  of  parameters  and  using  the  identical  approach 
mentioned  above,  excellent  FPC  performance  was  obtained. 
The  comparison  between  the  PI  and  FPC  control  responses  is 
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Figure  5.11  Comparison  Of  Simulation  Responses  For  PI  And 

Multiple  Pass  FPC  Control  Of  Bottom  Composition 
For  a  20%  Step  Decrease  In  Feed  Flow  Rate 
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shown  in  Figure  5.12.  The  ISE  value  has  decreased  from 
1.575  to  0.835. 

5.3.3  SET  POINT  CONTROL 

The  control  responses  for  a  2%  step  increase  of  set 
point  under  FPC  operation  using  the  parameters  estimated  in 
Run  B12  has  not  performed  well,  but  the  ISE  result  is  still 
comparable  to  that  achieved  using  a  well -tuned  PI 
controller.  It  is  not  surprising  that  the  FPC  control 
response  is  not  optimum  since  the  parameters  used  were  those 
estimated  from  load  disturbances.  From  the  responses 
plotted  in  Figure  5.13,  it  can  be  seen  that  both  the  input 
and  output  variables  for  the  FPC  run  oscillate  throughout 
the  entire  duration  of  the  run. 

To  improve  the  control,  it  was  decided  to  adapt  the 
parameters  obtained  in  Run  B12.  Allowing  the  parameters  to 
adapt  for  a  set  point  change  has  resulted  in  an  improvement 
of  control  performance  as  shown  by  the  lower  ISE  value.  The 
improved  control  performance  is  obvious  from  the  responses 
plotted  in  Figure  5.14,  although  a  slight  oscillation  in  the 
composition  response  still  remains.  However,  by  using  the 
final  parameters  obtained  in  this  run  and  adopting  the  FPC 
approach,  the  ISE  is  further  reduced  and  the  slight 
oscillation  is  eliminated  as  can  be  seen  from  the  responses 
in  F igure  5.15. 

For  the  case  of  a  2%  step  decrease  in  the  set  point, 
using  the  FPC  approach  with  the  parameters  obtained  from  a 
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Figure  5.12  Comparison  Of  Simulation  Responses  For  PI  And 

Multiple  Pass  FPC  Control  Of  Bottom  Composition 
For  a  20%  Step  Increase  In  Feed  Flow  Rate 
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Figure  5.13  Comparison  Of  Simulation  Responses  For  PI  And 

Multiple  Pass  FPC  Control  Of  Bo t tom  Compos i t ion 
For  a  2%  Step  Increase  In  Set  Point 
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p,jQjjP0  ^  ^ f\  S i mu  1  a t ed  Response  Of  STR  Control  Of  Bottom 

product  Composition  For  a  2%  Step  Increase  In 
Set  Point,  Starting  With  Initial  Parameters 
Values  As  Specified  In  Table  5.4 
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Figure  5.15  Simulated  Response  For  FPC  Control  Of  Bottom 

Product  Composition  For  a  2%  Step  Increase  In 
Set  Point,  With  Parameters  Obtained  From  Run  B20 
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-20%  feed  flow  rate  disturbance,  the  control  performance  is 
better  than  that  achieved  using  PI  control.  The  controlled 
responses  for  both  PI  and  FPC  are  plotted  in  Figure  5.16. 

By  all  owing  the  parameters  to  adapt  as  done  in  the  case  of 
the  increase  in  set  point,  further  improvement  in  control 
performance  and  reduction  in  ISE  value  were  accomplished. 
The  final  FPC  performance  is  shown  in  Figure  5.17.  A 
summary  of  FPC  and  STR  tests  performed  for  set  point 
disturbances  is  provided  in  Table  5.4. 
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Figure  5.16  Comparison  Of  Simulation  Responses  For  PI  And 

Multiple  Pass  FPC  Control  Of  Bottom  Product 
Composition  For  a  2%  Step  Decrease  In  Set  Point 
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Figure  5.17  Simulated  Response  For  FPC  Control  Of  Bottom 

Product  Composition  For  a  2%  Step  Decrease  In 
Set  Point,  With  Parameters  Obtained  From  Run  B23 
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5.4  Summary  Of  FPC  And  STR  Runs  For  Bottom 
Product  Composition  Control  For  Set 
Point  Disturbances 
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5.3.4  Effect  Of  Time  Delay  Tests: 

For  bottom  composition  control,  the  large  time  delay  is 
the  main  reason  that  the  control  performance  is  so  poor. 

The  effect  of  time  delay  on  the  self-tuning  regulator  has 
been  discussed  theoretically  by  Astrom  (5).  The  variance  of 
the  output  variable  is  expected  to  increase  when  there  is 
added  time  delay  in  the  system. 

To  investigate  how  system  time  delay  affected  control 
system  performance,  the  inserted  time  delay  was  changed  from 
16  minutes  to  4  minutes.  The  multiple  pass  method  was  used. 
The  parameter  estimation  runs  were  performed  using  the  same 
procedure  as  outlined  in  Section  5.4.1.  Instead  of  having 
to  estimate  7  parameters,  due  to  the  reduced  time  delay  that 
resulted  in  a  change  of  the  K  value  from  4  to  1 ,  only  4 
parameters  are  now  required.  The  runs  showing  the  effect  of 
the  reduced  time  delay  on  the  ISE  and  the  estimated 
parameter  values  are  summarized  in  Table  5.5. 
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Table  5.5 


Summary  Of  Reduced  Time  Delay 
Runs  For  Bottom  Product  Composition 
Control  Under  FPC  Operation 
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By  reducing  the  time  delay  from  16  minutes  to  4 
minutes,  the  control  performance  for  the  feed  flow 
disturbances  is  excellent.  The  ISE  values  decreased  from 
0.792  to  0.115  for  a  negative  feed  flow  rate  disturbance, 
and  from  0.835  to  0.108  for  the  positive  feed  flow  rate 
disturbance.  In  the  case  of  a  step  increase  in  set  point, 
the  ISE  values  decreased  from  0.083  to  0.026.  However,  for 
the  set  point  decrease  of  2%  test,  use  of  the  FPC  parameters 
resulted  in  a  performance  that  caused  the  ISE  value  to 
increase  from  0.073  to  0.1215.  The  fact  that  these  FPC 
parameters  were  adapted  for  feed  flow  disturbances  may 
account  for  the  increase  in  ISE  for  the  set  point  run. 

Reducing  the  time  delay  of  the  experimental  equipment 
is  possible  by  shortening  the  bottom  measurement  system 
recycle  line  and  improving  the  G.C.  filtration  system.  By 
model  formulation,  the  reduction  of  K  value  required  in  the 
STR  can  also  be  accomplished  by  an  increase  of  system 
sampling  interval.  The  results  from  the  tests  performed  for 
different  sampling  intervals  are  given  in  Appendix  H. 


5.4  Discussion 

The  simulation  results  for  proportional  plus  integral 
control  have  demonstrated  the  difficulty  in  satisfactorily 
controlling  bottom  product  composition.  The  long  time  delay 
and  time  constant  clearly  dominate  the  output  responses. 

The  use  of  only  a  feedback  control  scheme  has  no  capability 
to  compensate  for  the  effects  of  the  long  time  delay  since 
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no  corrective  action  is  possible  until  the  output  variable 
deviates  from  the  set  point.  Consequently,  to  avoid 
excessive  control,  the  controller  setting  has  to  be  eased 
off  and  the  control  behavior  becomes  sluggish. 

When  zero  initial  parameter  estimates  were  used,  the 
self- tuning  regulator  did  not  perform  well  either.  By  . 
adopting  the  multiple  pass  method,  the  necessary  adaption 
for  the  parameter  estimates  is  obtained.  The  fact  that  the 
model  has  a  correct  time  delay  parameter  defined  provides  an 
inherent  predictor  to  cope  with  the  time  delay  problem.  As 
a  result,  the  multiple  pass  approach  yielded  superior 
control  when  compared  with  well  tuned  PI  control  runs.  The 
acceptable  control  performance  of  the  multiple  pass  approach 
for  set  point  disturbances  using  parameters  that  were  tuned 
for  load  disturbances  has  demonstrated  the  robustness  of  the 
control  scheme.  If  the  parameters  are  allowed  to  adapt, 
even  better  performance  can  be  obtained.  Furthermore,  the 
control  performance  can  be  greatly  improved  by  reducing  the 
system  time  delay,  providing  the  appropriate  set  of  system 
parameters  are  used  in  the  multiple  pass  method. 
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6.  CONCLUSIONS  AND  RECOMMENDATIONS 
Control  of  top  product  composition  using  the  self-tuning 
regulator  has  shown  that  the  top  product  composition  control 
system  is  fast  and  stable.  This  is  due  to  the  fact  that  the 
loop  contains  only  a  short  time  delay  and  is  relatively 
noise  free.  A  PI  controller,  once  well -tuned,  also  provided 
excellent  performance  in  both  simulation  and  experimental 
tests.  The  PI  controller  outperformed  the  STR,  using  zero 
initial  parameters.  In  contrast,  the  bottom  product 
composition  loop  contains  a  long  time  delay  so  the  PI 
controller  was  not  able  to  provide  satisfactory  control. 

The  STR,  with  the  advantage  of  an  inherent  predictor  in  the 
model  formulation,  managed  to  achieve  a  much  improved 
control  performance  once  reasonable  estimates  were  obtained 
with  the  multiple  pass  method. 

The  basic  self-tuning  regulator  provides  satisfactory 
control  behavior  for  step  changes  in  feed  flow  rate  and  set 
point  even  though  the  regulator  is  not  formulated  for  such 
changes.  The  regulator  was  able  to  adapt  the  parameters  for 
the  disturbance  to  perform  adequately.  Tests  using  the 
fixed  parameter  controller,  despite  the  fact  that  the 
parameters  were  identified  for  a  different  disturbance, 
showed  that  the  controller  was  robust  enough  to  provide  good 
overa 1 1  control . 

On  the  basis  of  the  results  from  the  series  of  tests 
conducted  on  the  effect  of  STR  design  parameters  on  control 
performance,  the  following  conclusions  can  be  stated! 
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a }  Seal ing  factor ,  ^  : 

The  choice  of  the  value  significantly  influences  the 
control  behavior.  A  decrease  of  the  value  has  the  effect 
of  increasing  the  controller  gain,  as  a  result  the  change  in 
controller  output  at  each  control  interval  becomes  larger. 
Consequently,  the  parameters  converge  at  a  faster  rate.  In 
addition,  the  values  of  the  final  converged  parameters 
change  with  different  values.  The  smaller  the  value 
used,  which  causes  an  increase  in  the  difference  between 
and  -1.0,  and  a  decrease  in  the  absolute  value  of  fo^;,  the 
1 arger  the  offset . 

b)  Initial  covariance  matrix,  £(0): 

The  use  of  a  high  P(0)  value  which  implies  the  lack  of 
confidence  in  the  initial  parameter  estimates,  allows  a 
larger  change  in  the  values  of  the  estimated  parameters  at 
each  control  interval.  Tests  using  zero  initial  parameter 
estimates  showed  that  use  of  a  higher  value  resulted  in 
improved  control  performance.  By  starting  with  a  high  P.(0) 
value,  showing  less  confidence  in  the  parameters,  there  is 
more  chance  to  converge  to  the  true  values.  Consequently, 
the  offset  was  found  to  decrease  with  an  increase  in  the 
P ( 0 )  value.  However,  too  high  a  P(0)  value  causes  a  very 

r?  s: 

large  fluctuation  in  the  estimated  parameters  at  the  time 
the  disturbance  is  introduced  and  can  lead  to  an  unstable 


behavior . 


■ 
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c)  The  number  of  parameters,  L: 

Limited  simulation  tests  performed  indicate  that  if  a 
higher  than  an  adequate  number  of  ^  parameters  is  used,  only 
marginal  improvement  may  be  realized. 

d)  Initial  parameters  estimates,  8(0): 

Using  zero  initial  parameters  estimates  and  a  high  P_(0) 
value  provides  satisfactory  top  product  composition  control. 
For  bottom  product  composition  control,  the  STR  has 
difficulty  in  identifying  the  system  parameters  when  zero 
initial  parameter  estimates  are  used.  This  is  due  to  the 
large  time  delay  relative  to  system  time  constant. 

Parameter  estimation  runs  were  required  to  provide  good 
initial  system  parameters  for  subsequent  satisfactory 
control  of  bottom  product  composition. 


Recommendations  for  future  studies  are  the  following: 

a)  The  inserted  time  delay  in  the  simulated  bottom 
composition  control  loop  reflects  the  existing  time  delay  in 
the  actual  distillation  column.  The  major  problem 
encountered  in  the  control  of  bottom  composition  was  the 
long  time  delay.  The  reduction  of  experimental  system  time 
delay  can  be  accomplished  by  redesigning  the  G.C.  sample 
circulation  system  and  an  improvement  and  size  reduction  of 
the  filtration  equipment.  The  G.C.  analysis  time  could 
liKely  be  improved  by  experimenting  with  column  packing 


•  *, 
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other  than  Poropak  Q. 

b)  Adding  control  weighting  to  the  regulator  as  proposed  by 
Clarke  and  Gawthrop  (15)  is  recommended  to  remove  the  bang 
bang  performance  that  is  evident  in  this  work. 

c)  Incremental  control  output  [ u ( t )  -  u ( t - 1 ) ]  in  the  5TR 
formulation  to  remove  offset  should  be  included. 


Nomenclature 


Alphabet ic 


A 

a 

B 

b 

C 

c 

E 

e 

F 

G 

H 

H 

K 

L 

n 

E 

t 

u 

V 

X 


Polynomial  in  Astrom  model 
Coefficient  in  polynomial  A 
Polynomial  in  Astrom  model 
Coefficient  in  polynomial  B 
Polynomial  in  Astrom  model 
Coefficient  in  polynomial  C 
Expectation  operator 
Normal  random  variable 
Polynomi a  1 
Pol ynomi a  1 
Vector 
Gain  vector 
Time  delay 

Number  of  parameters 
Order  of  Model 
Covariance  matrix 
sampling  instant 
Control  vector 
Loss  function 
Output  vector 


Greek 

o(L 


Coefficient  of  predictive  model 


, 


j  r  ,  s<  sot  ^HsvoO 


. 
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ent  of  predictive  model 
factor 

verage  process  of  driving  noise  e 
r  vector 

ial  forgetting  factor 
ion  vector 


ement 

onal  controller  gain  value 
me 

limit 

d)  Superscripts 

T  Matrix  transpose 

-1  Matrix  inverse 

^  Estimated  value 

e)  Abbreviation 

DDC  Direct  digital  control 

GC  Gas  Chromatograph 

PI  Proportional  plus  integral 

STR  Self-tuning  regulator 


p* 

Coef f ici 

b 

Seal ing 

e 

Moving  a 

a 

Par amete 

\ 

Constant 

/* 

Exponent 

± 

Informat 

Subscripts 

• 

i ( th)  el 

kp 

Proport i 

Tr 

Reset  ti 

UL 

Control 

not osl  ffif  ii  .ent  qx3 
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APPENDIX  A:  RECURSIVE  LEAST  SQUARES  ESTIMATION 


In  the  self -tuning  regulator,  a  recursive  least 
squares  algorithm  is  used  to  update  the  parameters 
at  every  sampling  instant  based  on  current  and  previous 
input  and  output  data.  The  equations  used  are  A.1, 

A. 2  and  A. 3. 


f(t)  =  f(t-1)  +  K  ( t )  |  y  ( t )  -  j?0u(t-k-1)  - 

^(t-k-ueU-i)] 


(A.1) 


K(t)  =  £(t-1)i(t-K-1)|1+4T(t-K-1)Plt-1)i(t-K-1)] 
£(  t )  -  K  ( t )  [  1+±T(t-k-1  )P(t-1 ) 


(A. 2) 


±(t-k-D)  KT(  1 ) } 


The  updating  gain  matrix  K(t)  in  Equation  A. 2  is 
calculated  first,  then  the  covariance  matrix  P(J^)  in 
Equation  A. 3,  followed  by  the  parameter  vector  0(t) 
i n  Equa t ion  A.1. 

The  recursive  least  squares  algorithm  is  derived  in 
the  following  manner: 

The  predictive  model  used  for  estimating  the  self-tuning 
regulator  parameters  on-line  is  given  by  Equation  A. 4: 


y  ( t )  +«f1y(t-k-1 )  +  cX^y  l  t-k-2 )  +  ...  tAft-k-n) 

=  (-Li  u(t-k-l)  +  (?, u ( t -k-2 )  +  ...+  fl  u(t-l-k-L)I+e(t) 

1  •“  (A. A) 


where  L  =  n+K-1 


Defining : 


±Tm  =  I  -y(t),  -y(t-1), 


•  •  •  • 


-y ( t -ri+ 1 ) , 


( A .  5 ) 


^  ( t )  -  I  DC  ^  •  (XCj » 


•  •  •  » 


£tn  '  P’l  ’  Pl  ^ 


(A. 6) 


■ 
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all  ow  s  Equation  A  .  A  to  be  w  r  i  1 1  e  n  as’ 


y(t)  =|?£>u(t-k-i)  +  'jfT( t -K- 1 )  e  4  £(t) 


(A. 7) 


Adopting  the  notation  used  by  Sage  £  Me  Isa  (7),  let 


Zt  =  y(t)  "  f20u(  t-k-  1 ) 

Ht=  iT(t-K-1) 

x.  =  e 

vt=  £(  t ) 

zc=  HtX  +  Vt 

The  least  squares  estimate  ofiS^X^is  as  follows: 

%=  it-i+  Mzt- 

!<t= 

£c=  £t-i  -  < 1  +  UtEt-i 


(A. 8) 

(A. 9) 
(A. 10) 
(A. 11) 


At  time  t,  are  available,  so  P_^-  is 

calculated  using  lPr  >  and  presently  a va i  1  ab  1  ej^  before  Jw- 
can  be  calculated?  L 


Another  form  of  Equation  A. 10  can  be  written  by 
substituting  Equation  A. 11  into  Equation  A. 10. 

5t=  [£t-l  -  £*-iHT<1  +  h^p^h/)-1 

=  P*-i  H?  +  r1 

=  P*.lUx  II  -  u  +  iUPt-iH?  ] 

=  £t-lH^[1  -  (H^P^jj/)  /  (1  +  id? )  1 

=  Lt-iui  [i  +  (A 

For  this  study,  Equation  A. 12  is  used  instead  of 
Equation  A. 10.  would  be  calculated  at  time  t  and  used 

at  time  T>1  foF  calculating  K^^to  obtain  X^^.. 


:s..«T  a*  *«  ai  »**»■«*•  «-»*upa  tssst  art' 
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The  chromatograph  analysis  of  samples  of  the  bottom 
product  provides  the  composition  measurement.  The  gas 
chromatogr aph  analysis  time  and  sample  system 
transportation  time  delay  plus  the  process  time  delay 
constitute  the  total  system  time  delay. 


1.  Gas  chromatograph  analysis  time: 

The  gas  chromatograph  cycle  time  limits  the 
sampling  period  to  no  shorter  than  256  seconds  due 
to  the  PoropaK  Q  column  packing  conditions  required 
for  a  reliable  analysis. 


2.  Transpor ta t ion  time  delay  calculation: 

The  product  circulation  system  arrangement 
contributes  a  transport  time  delay  of  up  to 
13  minutes.  The  system  configuration  that  results 
in  this  time  delay  is  shown  schematically  in  Figure  B1. 
The  magnitude  of  the  time  delay  was  calculated 
as  follows: 


Section  1:  Distance  between  reboiler  bottom  and 

pump  (cf  Figure  B1)  =  210.00  cm 

Size  of  tubing  (copper  type  L) 

0 . D .  =  0.95  cm 

I . D .  =  0.79  cm 

•  Thickness  =  0.08  cm 

3 

Based  on  a  volumetric  flow  rate  of  0.83  cm  /sec., 
the  transport  time  delay,  excluding  the  effect 
of  the  heat  exchanger  is: 

=  (distance)  *  ( cross-sect iona 1  area) 

(volumetric  flowrate) 

=  123.00  seconds 


T  ime 


■ 
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Section  2:  Distance  between  pump  and  second 

filter  =  486.41  cm 

Size  of  tubing  (copper  type  L) 

0 .  D  .  -  0.64  cm 
1 . D .  =  0.48  cm 
Thickness  =  0.08  cm 

Velocity  =  4.10  cm/sec 

Time  to  travel  specified  distance  =  118.64  seconds 


Section  3:  Distance  between  second  filter  and 

gas  chromatograph 

Size  of  tubing  (plastic) 

I . D .  =  0 . 25  cm 

Velocity  =  14.60  cm/sec 

Time  to  travel  specified  distance  =  6.52  seconds 


Total  transportation  time  for  all  three  sections 
=  248.16  seconds 


Contribution  of  heat  exchanger  to  transportation 
time  delay: 

The  heat  exchanger  is  located  between  the  reboiler 
and  pump  as  shown  in  Figure  B1. 

Exchanger  manufacturer: 

Heat-X-Changer  Co.  Inc.,  Brewster,  New  York. 

Type:  Single  pass  finned  heat  exchanger. 

Speci f icat ion: 

0 .  D .  =  6 . 35  cm 

Length  =  30.96  cm 
Copper  pipe  fitting  O.D.  =  3.48  cm 
Copper  pipe  fitting  I.D.  =  3.18  cm 

Volume  of  fluid  in  heat  exchanger  to  be 
displaced  before  a  fresh  sample  reaches  the 
gas  chromatograph: 

Cross-sectional  area  of  heat  exchanger  =  7.92  cm 
Volume  of  heat  exchanger  =  7.92  cm  x  30.96  cm 

3  =  245.20  cm3 

Flowrate  =  0.83  cm  /sec 

Time  =  295.42  seconds 
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Contribution  of  filters  to  transpor tat  ion  time  delay: 

There  are  two  filters  in  series: 

First  filter  volume  =  not  significant 
Second  filter  volume  =  200.00  cm3 

Flowrate  =  0.83  cm3 /sec 

Time  =  240.96  seconds 


Total  time  delay  due  to  t r anspor t at i on : 
(248.16  +  295.42  +  240.96)  =  784.54  seconds 


3.  Process  time  delay  measurement: 


Experimental  time  delay  tests  with  samples 
taken  at  the  bottom  of  the  reboiler  and  downstream 
of  the  gas  chroma togr aph  heat  exchanger  provide  a 
reasonable  estimate  of  the  process  time  delay. 
Seven  samples,  at  256  second  intervals,  were  taken 
with  the  steam  flow  rate  decreased  by  5%  at  128 
seconds  of  elapsed  time.  The  results  are  given  in 
Figure  B2.  The  composition  of  the  samples  taken 
at  the  bottom  of  the  reboi ler  began  to  increase 
after  one  sample  interval  and  those  at  the  heat 
exchanger  did  not  change  for  two  sample  intervals 
due  to  extra  transport  time  delay.  Consequently, 
the  process  time  delay  was  estimated  as  one  sample 
interval,  256  seconds. 


As  can  be  seen,  the  total  time  delay  appears 
to  be  about  21  minutes,  so  with  a  sampling  interval 
of  256  seconds,  the  resulting  k  value  is  5.  This 
value  for  the  time  delay  is  consistent  with  that 
established  by  calculation. 
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Figure  31  •  Schematic  diagram  of  gas  chromatograph  sampling  system 
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Figure  B2  Dynamic  response  of  bottom  composition  to  a  steam  plow 

RATE  STEP  CHANGE  OF  -5%  WITH  SAMPLES  TAKEN  AT  REBOILER 
BOTTOM  AND  HEAT  EXCHANGER  OUTLET 
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Appendix  C:  Experimental  Control  Program  Listing 


1 

2 

3 

4 

5 

6 
7 
B 
9 

10 

11 
12 

13 

14 

15 

1 6 
17 
IB 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 


//  FOR  KKTOI 
• J  DCS (DISK) 

*  )DCS(PAPERTAPE  ) 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


LT 


SIR  CONTROL  PROGRAM 


ALGORITHMS  AVAILABLE 
1 .STR-U  METHOD 
2 . STR-DELU 
3. PI  CONTROL 


DISTURBANCE  USED 

1 .DETERMINISTIC 
2 . STOCHASTIC 


REAL  LTS.MU. LDUMT( 10).LDUM3( 10) . KPT . KI T . KPB , K IB 
DIMENSION  PT ( 8 , 8 ) , PB ( 8 , 8  )  , X  1  ( 8 ) , X2 ( 8 ) . XDUMT ( 10) ,XDUMB( 10) 
DATA  KPT .KIT, KPB. K1B/0. 0.0. 0,0. 0,0.0/ 

DATA  MT , M3 , KT , KB/2 ,2,1, O/ 

DATA  BOB , BOT , RT , RE , KST R/O .0.0. 0,0. 0.0. 0,0/ 

DATA  ACCMT, ACCM3. REST. RE SB/O. 0.0. 0.0. 0.0.0/ 

.DEFINE  FILE  50(  4 , 320.  U .  1  F I  LE  ) 

CALL  GETTV(LT) 

LUNV=LT 

IN  THIS  PROGRAM  IS  LUN  NUM3ER  BUT  IN  RKT02  IT  IS  FLOW 


403 


INITIALIZE  INPUT 
DO  403  1=1.10 
LDUMT (I ) =0 . O 
LDUM3 ( I )=0.0 
XDUM3( I ) =0.0 
XDUMT ( 1 )=0.0 
I BT  C=0 


AND  OUTPUT  DUMMY  VECTORS 


36 

1 PRT  =0 

37 

KK=0 

38 

C 

39 

C 

CONTROL  LOO 

40 

c 

FEED, TOP  CO 

41 

JFD= 1537 

42 

I XT  = 1540 

43 

1 LT  = 1 538 

44 

IXB= 1543 

45 

ISM« 1539 

46 

c 

47 

c 

ITEM  KEYS  1= 

48 

11  =  1 

49 

12  =  2 

50 

13  =  3 

51 

c 

52 

c 

SET  INITIAL 

•  53 

c 

SET  INITIAL 

54 

c 

INITIAL  PAR/ 

55 

DO  20  U= 1 , 8 

56 

DO  20  1*1.8 

57 

PB ( I . J)=0.0 

58 

20 

PT ( I ,d)=0-0 

59 

DO  2  1  I  =  1  .B 

60 

xi(i)=o.o 

I  .D. 

.REFLUX, 


BOTTOM  COMP.. STEAM 


1  =  ME  AS .  2  =  SET  PT  3  =  DUTPUT 


COVARIANCE 
COVARIANCE 
METERS  SET 


MATRIX  FOR 
MATRIX  FOR 
TO  ZERO 


TDP:  lOO. 
BOTTOM  lO. 
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61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

1O0 

101 

102 

103 

104 

105 

106 

107 

103 

109 

1  lO 

111 

112 

113 

114 

1 15 

1 16 

1 17 

1  18 

1  19 

120 


136 


X2(l )=0.0 
PT(  1 , 1 )= 100.0 
21  PB(1.!)*10.0 

C 

KU  = 1.0 
MVR«=0 
KPAR=4 
C 

WRIT  E ( LT , 1 1  ) 

11  F  ORMAT (  '  RUN  TIME  (  I  MTS) . PRINT  CNTRL(KPRT). BOTMS  SAMPL1 NG( 1 SMP )  '  ) 

CALL  FF INP(LT .3,0. IMT 5 .O, KPTR ,0. 1 SMP. IER) 

C 

C  GET  CURRENT  FLOWS  SET  PTS  AND  COMP.  MEAS. 

IT=  1 

CALL  GTVLU(IXT. I1.XTS.IER.IT) 

CALL  GTVLU( 1XB. I  1 .XBS, IER.  IT) 

IT  =  4 

CALL  GTVLU( IFD. I2.FD.  IER.  IT) 

CALL  GTVLU(ISM. I2.SFS. IER.  IT) 

CALL  GTVLUdLT.I2.LTS. IER, IT) 

WRIT E ( LT .13)  XTS. LTS.XBS, SFS, FD 

13  FORMAT ( / , ' T  C0Mk='F6.3.'  RFLX='F6.3,'  B  C0MP='F6.3,'  STM='F6.3, 

1'  F EED=  '  F6 . 3 ) 

C 

C  GET  CURRENT  OUTPUTS  IN  %  FROM  REFLUX  AND  STEAM 

CALL  GTVLU(ILT. I3.LTS.IER. 1) 

CALL  GTVLUdSM.  13.  SFS.  IER,  1) 

C 

C  SET  DISTURBANCE  TYPE 

WRITE  (LT.15) 

15  FORMAT ( 'KTYP--T  COMP  SP . 2  B  COMP  SP . 3  FEED  CHANGE  AND  MVR - -  O  STR 
1  1  MVR '  ) 

CALL  FF1NP(LT.2.0.KTYP.0.MVR. IER) 

C 

c 

C  UDEL=1  STR-U  JDEL=2  STR-DELU 

WR I T  E ( LT , 27 ) 

27  FORMAT( ' ENTER  UDEL  1 -STR-U . 2  =  STR -DELU  UDIS  1  =DET .  2  =  ST0C./) 

CALL  FFINP(LT.2.0.UDEL .0. JDIS. IER) 

C 

c 

c 

C  GET  TIME  AND  MAGNITUDE  OF  DISTURBANCE 

WR I T  E ( LT , 10) 

10  FORMAT(  'DISTURBANCE  F ROM ( KDI S )  AND  VALUE ( DI ST )' ) 

CALL  FFINP ( LT , 2 . O. KDI S , 1.D1ST.IER) 

WRITE  (LT.16) 

16  FORMAT ( ' KCONT  = 1 -T  STR  B  0PN.2-T  PI  B  0PN.3-T  OPN  B  STR.4-T  OPN  B  P 

II') 

WRITE(LT .26) 

26  F  ORMAT (  '  5-T  STR  B  PI.6-B0TH  PI.7-B0TH  STR.8-T  PI  B  STR') 

CALL  FFINP(LT. 1 ,0. KCONT, IER) 

GO  TO  (224. 1  13.225,  114.223, 111.222,  112) .KCONT 
113  WR ITE(LT,22) 

22  FORMAT( 'SPECIFY  KPT, KIT') 

CALL  FFINP(LT.2. 1 .KPT, 1 .KIT, IER)  . 

WRITE ( LT , 36  )  KPT. KIT 
CALL  MMANL(ILT) 

GO  TO  333 
224  WR I T  E ( LT , 23 ) 


- 


• 
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141 
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143 
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172 
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176 

177 
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23  FORMAT ( 'SPECIFY  BOT.RT,  JDENT  STOP  K S T R ' ) 

GALL  FF 1NP(LT .3. 1 .BOT. 1 .RT .O.KSTR. IER) 

VR I T  E ( L  T . 99 ) 

99  f0RMAT(  'SPECIFY  X1(I),  1=1,5') 

CALL  FFIUP(LT.5. 1 . X 1 ( 1). 1,X1(2).  1,X1(3). 1 . X 1 (4 ) . 1 , X  1  ( 5  )  . I ER) 

WRITE(LT .39)  MT.KT.BOT 

WRITE(LT.53) 

53  FORMAT( 'SPECIFY  X1(I). 1*6.8) 

CALL  F  F INP( LT.3, 1  ,  X  1  (  6  )  . 1 . X 1 ( 7 ) ,  1.X1(8),IER) 

WRI  TE ( LT . 98  )  ( X 1 ( I ) . I  =  1 . 5 ) 

WRIT  E(LT,92)(X1(I )  . I =6 . 8  ) 

92  FORMAT (5X. 'X1(G)=' .F7.4, 'X1(7)=' .F7.4. 'X1(8)=' .F7.4) 

CALL  MMANL(ILT) 

GO  TO  333 
225  VR I T  E ( LT , 24 ) 

24  FORMAT( ' SPECIFY  BOB  , RB . KB . I DENT  STOP  KSTR ' ) 

CALL  FFINP(LT. 4.1. BOB. 1.RB.0. KB. O.KSTR, IER) 

VRITE(LT,97) 

97  FDRMAT( 'SPECIFY  X2 (I). 1=1,5') 

CALL  F F I NP ( LT , 5 . 1 , X2 ( 1 ) . 1 , X2 ( 2 ) . 1 . X2 ( 3 ) . 1 . X2 ( 4 ) . 1 . X2 ( 5 ) . I ER ) 

WRI  T  E  (  LT  ,  54  ) 

54  FORMAT( ' SPECIFY  X2(I),I=6.8) 

CALL  FFINP(LT«3,1.X2(6). 1 . X2 ( 7 ) . 1 . X2 ( 8 ) . I ER ) 

WRITE (LT .41  )  MB, KB. BOB 
WRITE (LT ,96)(X2(I ).I=1 .5) 

WRITE (LT,93)(X2(I).I=6.8) 

93  FORMAT (5X.  ' X2 ( 6  )  = ' . F7 . 4 .  ' X2 ( 7 )  =  ' . F7 . 4 .  ' X2 ( 8 ) = ' . F7 . 4 ) 

CALL  MMANL(ISM) 

CALL  MMANL(IXT) 

GO  TO  333 
114  WRITE(LT,25) 

25  FORMAT ( 'SPECIFY  KPB.KIB') 

CALL  F F 1 NP ( LT ,2 , 1.KPB, 1.KIB.IER) 

WWI T E ( LT ,48)  KPB.KIB 
CALL  MMANL(ISM) 

CALL  MMANL(IXT) 

GO  TO  333 
111  WRITE (LT. 17) 

17  FORMAT( 'SPECIFY  KPT . KI T . KPB . K IB ' ) 

CALL  FFINP(LT ,4. 1 .KPT . 1 .KIT, 1 .KPB. 1 ,KIB. IER) 

WRITE (LUNW. 40) 

WRIT  E ( LUNW ,42)  KPT . K I T . KPB . KI B 
CALL  MMANL(ILT) 

CALL  MMANL(ISM) 

GO  TO  333 
222  CONTINUE 

WRITE  ( LT .14) 

14  F ORMAT (  'BETA  O  ( BOT , BOB ) ,  NOISE  COVAR ( RT , RB ) . KB . 1  DENT  STP(KSTR)') 

CALL  F  F I NP ( LT , 6 , 1 . BOT .  1 .BOB.  1 . RT . 1 . RB . O . KB . O. KSTR . I ER ) 

WRITE(LT .99) 

CALL  FFINP(LT.5.1.X1(1).1,X1(2).1.X1(3),1.X1(4).1,X1(5).IER) 

WR 1  T  E ( LT . 53  * 

CALL  F  F 1 NP (LT.3, 1,X1(6), 1,X1(7),1.X1(8),IER) 

WRI TE ( LT , 97  ) 

CALL  FF 1 NP ( LT.5. 1 ,X2( 1). 1.X2( 2). 1.X2( 3), 1.X2 (4). 1,X2( 5). IER) 

WRI TE ( LT . 54  ) 

CALL  FFINP(LT,3. 1.X2(6). 1 . X2 ( 7 ) . 1 . X2 ( 8 ) . I ER ) 

WRITE  (LUNW, 35) 

WR 1 T  E ( LT . 37  )  MT . KT . MB . KB . BOT . BOB 
WR I T  E ( LT  , 98  )  ( X  1  ( I ) . I  *  1 . 5 ) 


# 


. 


- 
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201 
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204 
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207 
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WR 1 T  E ( L  T , 9  2 ) ( X  1  (  I  )  ,  I  c  6 . 8  ) 

WRITE(LT .96)  (X2( 1  )  .  1  -  1.5) 

WR1TE(LT.93)(X2(1 ),I=6,8) 

CALL  MMANL(ILT) 

CALL  MMANL(ISM) 

GO  TO  333 

1 12  WRITE  (LT. 18) 

18  FURMAT( 'SPECIFY  KPT. KIT. AND  BOB . RB . KB ' ) 

CALL  FFINP(LT,5,  1.KPT,  1.KIT.1  .BOB. 1 . RB . O . KB , I E R ) 

WRITE(LT,97) 

CALL  FFINP(LT .5,  1 ,X2(  1  ).  1  ,X2(2), 1 , X  2 ( 3 ) ,  1 ,X2(4). 1 ,X2(5), IER) 

WRITE (LT .54) 

CALL  F  F I NP ( LT . 3 .  1  . X2 ( 6  )  . 1 , X2 ( 7 ) . 1 . X2  (  8 )  .  I ER ) 

WRITE  ( LT .44) 

WR I T E ( LT .45)  KPT . KI T , MB . KB . BOB 

WR I T  E ( LT , 96 ) ( X2 ( I ) . I  =  1 , 5 ) 

WRITE(LT.93)(X2(I),1=6.8) 

CALL  MMANL(ILT) 

CALL  MMANL(ISM) 

GO  TO  333 
223  WR ITE(LT, 19) 

19  FORMAT( ' SPECIFY  BOT.RT.  AND  KPB, KIB') 

CALL  F  F I NP ( LT .4,1 . BOT .  1  .  RT . 1 ,KPB.  1  .KIB. IER) 

WR I T  E ( LT . 99 ) 

CALL  FF1NP(LT,5, 1 , X 1 (  1).  1.X1(2).1,X1(3).1.X1(4),  1,X1(5),1ER) 
WRITE(LT.53) 

CALL  FFINP(LT,3,  1 , X 1 ( 6 ) ,  1,X1(7).  1,X1(8).IER) 

WRITE  (LT.46) 

WR ITE(LT,47  )  MT . KT . BOT . KPB  .  KIB 
WR I T  E ( LT , 98 )  ( X  1  ( I  )  . I  =  1 . 5 ) 

WRITE (LT .92) (XI ( 1 ) . 1=6.8) 

CALL  MMANL(ILT) 

CALL  MMANL(ISM) 

333  CONTINUE 

CALL  T 1  ME ( I  HR , 1 MT , I  SEC ) 

WR I TE ( LT , 38 )  IHR.IMT.1SEC 

1CYC=640 

REFT  =  XT  S 

REFB  =  XBS 

WRITE (50'  1 )  KPT R. REFT .REFB. 1 PRT , ACCMT , ACCMB , MU , KPAR . KSTR . MVR , KK 
1  .  JDEL . UDI S.MT ,KT .MB .KB . RT . RB. BOB, BOT  .  F  D 
WR I T  E ( 50 ' 2 )  LDUMT . LDUMB . XDUMT , XDUMB . X 1 . X2 . LT S , SFS 
WRITE (50' 3)  PT.PB 

WRITE (50' 4  )  KCONT , 1 MT S . LUNW . KOI S . DI  ST  .  KT YP . I  BTC. I SMP . XTS . XBS . 

1LTS. SFS. KPT. KPB. KIT. KIB. RE ST .RE SB 
I LEV=B 
IBIT=2 
I TMR=7 

C  TO  PUT  THE  TIMER  IN  RESTART  TABLE  ****•**•»•  ****** 

C 

CALL  SCPTB ( 1 , ITMR, 1  LEV, 1BIT, 1 . 1CYC.IE) 

CALL  CYCLE( ILEV, IBIT. ITMR. 1CYC) 

35  FORMAT  ( /  ,  lOX ,  ' *  *  STR  ON  BOTH  COMPOSITIONS**') 

36  FORMAT ( /  .  10X  ,  ' * • P 1  ON  TOP  AND  BOTTOM  COMP  OPEN* * ' . / . 5X ,  ' KPT= ' F8 . 4 . 
1 '  KI T  = ' ,F8 .4) 

37  FORMAT ( / , 5X ,  '  MT = '  . 1 2  .  '  KT=',I2.'  MB='.I2.'  KB='I2,'  BOT=',r8.4. 

1'  BOB  =  '  F  8 . 4  ) 

38  FORMAT(/.5X. 'CONTROL  ACTION  INITIATED  AT'. 313./) 

39  F  ORMAT ( /  .  lOX . ' * • STR  ON  TOP  COPMP . .  BOT  M  OPEN  LOOP* •'.//. 5X .' MT , 
112,'  KT  *',12.'  BOT  = ' , F  8 . 4 ) 


' 


- 
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24  1 
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247 
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252 

253 

254 

255 

256 

257 

258 

259 


40  F0RMAT(/, 10X. '  **PI  ON  BOTH  COMPOSITIONS**') 

41  F  ORMAT ( / ,  lOX  .  ' *  *  T  OP  OPEN  LOOP,  STR  ON  BUT  M  COMP .**'.//. 5X MB ®', I  2 
1. '  KB® ' ,12.  '  BOB®' , F  8 . 4 ) 

42  FORMAT (/ , 5X , 'KPT  = ' F8 . 4 , '  KIT»',F8.4.'  KPB='.F8.4.'  KIB=',F8.4) 

44  FORMAT (/ , 10X , ' * "PI  ON  TOP  AND  STR  ON  BOT  COMPS.**') 

45  F  ORMAT (/ ,5X,  ' KPT  =  ' FB . 4 ,  '  KIT=',F8.4.'  M8=',I2.'  KB®'. 12, 

1'  BOB  = ' , F8 . 4  ) 

46  FORMAT(/. 10X , ' * *STR  ON  TOP  AND  PI  ON  BDT  COMPS.**') 

47  FORMAT (/ , 5X ,  ' MT = ' . 1 2 .  ' KT® ' . 1 2 .  'BOT® ' ,FB. 4. ' KPB® '. FB. 4.  'KIB*'.F8. 4) 

48  F  ORMAT ( / .  10X  ,  ' *  *  T  OP  OPEN  LOOP.  PI  ON  BOT  M  COMP .**',//. 5X .' KPB ® ' 
1.F8.4.  '  KIB®' , F  8 . 4 ) 

98  FORMAT (/ , 5X ,  'X1(1)=',F7.4,'X1(2)='.F7.4.'Xl(3)=',F7.4.'X1(4)*'t 
1F7.4, *X1(5)  =  '  ,F7.4) 

96  F  ORMAT ( / , 5X ,  'X2(l)='.F7.4.'X2(2)='.F7.4.'X2(3)=',F7.4.'X2(4)='. 

1F7.4,  ' X2 ( 5 ) * ' , F7 . 4 ) 

CALL  EXIT 
END 

•DELET  A  RKTOI  »**»* 

®  STORECI I  A  1  RKTOI  RKTOI 


260 

•FILES(BO.DISTO.O) 

261 

•CCEND 

262 

// 

END 

263 

264 

// 

FOR  RKT02 

265 

*  I OCS ( D I SK ) 

266 

*  I OCS( PAPERTAPE ) 

267 

REAL  LTS . LT .KPT , KIT . KPB 

268 

DIMENSION  PT ( 8 , 8 ) , PB ( 8 , 

269 

DEFINE  FILE  50(4,320.0. 

270 

DEFINE  FILE  51(240, 46. U 

271 

IFD=  1537 

272 

I LT  =  1 538 

273 

ISM® 1539 

274 

I XT  *  1540 

275 

IXB= 1543 

276 

IBM= 1544 

277 

11  =  1 

278 

12  =  2 

279 

13  =  3 

280 

READ  (  50 '  1  )  KPTR  ,  REFT  .  Ri 

28  1 

1 , ODE  L . JD I S , MT . KT . MB . KB .  1 

282 

READ  (50' 2)  LDUMT.LDUMB 

283 

READ  ( 50 ' 3 )  PT , PB 

284 

READ  ( 50 ' 4 )  KCONT.IMTS.I 

285 

1LTS.SFS.KPT .KPB. KIT .KIB 

286 

KK=KK+ 1 

287 

1 PRT  =  I PRT  + 1 

288 

1 BT  C= 1BTC+ 1 

289 

IF( JDIS- 1 )23, 23.25 

230 

23 

IF(KK-KDIS)  333.14,333 

291 

14 

CONTINUE 

292 

GO  TO  (44,22.55),  KTYP 

293 

44 

XTS=DI ST 

294 

GO  TO  333 

295 

22 

XBS=DIST 

296 

GO  TO  333 

297 

55 

CALL  PTVLU( 1FD. I2.DIST,  ] 

298 

333 

CONTINUE 

299 

GO  TO  17 

300 

25  IF(KK-30) 17.27.27 

) 


.  r 


' 


- 


- 


' 
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301 

27 

CAU  GAUSS( 333.DI ST . 18. .SDIST) 

302 

CALL  PTVLU(IFD. 12. SDIST. IER) 

303 

17 

CONTINUE 

304 

C 

GET  THE  DATA  FROM  PROCESS  MODEL 

305 

C 

306 

17-1 

307 

CALL  GTVLUdXT.I1,  XT  .IER. IT) 

308 

CALL  GTVLU(IXB. I 1 .XB. IER. IT) 

309 

CALL  GTVLUdBM.11.  EM. IER. IT) 

310 

I  T*=  4 

311 

CALL  GTVLUflFD. I 1 . FD. IER. IT) 

312 

• 

IF(KK-IMTS)  11,12.12 

313 

11 

CONTINUE 

314 

CALL  GTVLUdLT. I  1. VLT.JER.IT) 

315 

CALL  GTVLU(ISM, 1 1 .SM. IER. IT) 

316 

GO  TO  (222.111.33.33.222.111.222.111).  KCONT 

317 

111 

CALL  GTVLU(ILT. 13, LT. IER.  1) 

318 

CALL  PID(KPT, KIT, REST. REFT. XTS, LTS, XT, CLT. ACCMT ) 

319 

CALL  PTVLU(ILT. I3.CLT. IER) 

320 

GO  TO  33 

321 

222 

CONTINUE 

322 

OLT=LT 

323 

CALL  GTVLUdLT. 13, LT. IER. 1) 

324 

XTN= ( XT -XT  S ) /R  E  FT 

325 

I F (UDEL- 1 )69,69,70 

326 

70 

LTN= ( LT-OLT )/LTS 

327 

OLT  =LT 

328 

GO  TO  71 

329 

70 

LTN=( LT-LTS )/LTS 

330 

71 

CALL  STR  (  MT . KT . BOT . RT . PT . X  1  . LDUMT . XDUMT . XTN . LTS . LTN . CLT , ACCMT . OLT ) 

331 

CALL  P7VLU( ILT, 13. CLT , I ER) 

332 

33 

CONTINUE 

333 

IF(IBTC-ISMP)  13,15.13  ~ 

334 

15 

I BTC=0 

335 

GOTO  (13.13,223.112,112.112.223.223),  KCONT 

336 

1  12 

CALL  GTVLU(ISM.I3.SF.IER. 1) 

337 

CALL  PI D( KPB , KIB , RE  SB .REFB.XBS.SFS.XB.CSF.  ACCMB ) 

338 

CALL  PTVLU(ISM. 13. CSF . IER) 

339 

GO  TO  16 

340 

223 

CONTINUE 

34  1 

OSF=SF 

342 

343 

344 

345 

346 

347 
34B 

349 

350 

351 

352 

353 

354 

355 

356 

357 

358 

359 

360 


67 


16 


12 

C 

C 


CALL  GTVLU( ISM. 13 . 5F . I ER, 1 ) 

XBN=-(BM-XBS)/REFB 
IF  (JDEL-1 >68,68,86 
86  SFN=(SF-OSF)/SFS 
OSF  =  SF 
GO  TO  67 

68  SFN=(SF-SFS)/SFS 

CALL  STR(MB .KB . BOB . RB.PB.X2  .  LOUMB .XDUMB . XBN , SFS . SFN . CSF . ACCMB.OSF ) 
CALL  PTVLU(1SM, 13. CSF. IER) 

CONTINUE 
GO  TO  13 
CALL  CANCL ( 07 ) 

TO  REMOVE  THE  TIMER  FROM  RESTART  TABLE  ••*••• 

CALL  SCPTR(0. 7 .8 . 2 . 1 .600. IE ) 

CALL  PTVLU(IFD. I2.FDS. IER) 

CALL  PTVLU(ILT . 13. LTS. IER) 

CALL  PTVLUdSM.  13.SFS.  IER) 

CALL  T 1  ME ( I HR , I MT , I  SEC ) 


t 


0 

3G1 

302 

363 
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368 
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371 
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375 

376 

377 

378 

379 

380 

381 

382 

383 

384 

385 

386 

387 

388 

389 

390 

391 

392 

393 

394 

395 

396 

397 

398 

399 

400 

401 

402 

403 

404 

405 

406 

407 

408 

409 

410 

4  1 1 

412 

4  13 

414 

415 

416 

417 

4  18 

419 

420 


v  ro 
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WRITE  (LUNW.3B)  1  MT . 1  SEC , XT .LT , XB . SF ,FD, ACCMT , ACCMB 
WRITE  (LUNV.66)  1HR.IMT.ISEC 
66  FORMAT (/, 'JOB  COMPLETED  AT  '.313) 

C  PUT  DATA  ACCUMULATION  LOOPS  NON-OPERABLE 

CALL  NONOP( 1584 ) 

CALL  NONOP ( 1585) 

CALL  NDNDP ( 1586 ) 

CALL  NONOP ( 1587) 

CALL  NONOP ( 1 5B8 ) 

CALL  N0N0P(1589) 

CALL  NONOP ( 1590) 

CALL  NONOP ( 1 59 1 ) 

CALL  NONOP ( 1592) 

CALL  NONOP (1667) 

GO  TO  999 
13  CONTINUE 

CALL  T I ME ( I HR . I MT . I  SEC ) 

1F( I PRT-KPTR )  412.413.412 
413  I PRT  =  0 

WRITE  (LUNW.38)  I MT . 1 SE C . XT , VLT . XB . SM , FD . ACCMT . ACCMB 
412  CONTINUE 

38  FORMAT (/. 21 3. ' XT= ' . F6 . 3 .  '  LT=  '  , F6 . 3 . ' XB= ' , F6 . 3 . 'SF=' ,F6.3, 

T ' FD= ' , F6 . 3 , 2F9 . 7  ) 

WR I T E ( 50 '  1 )  KPTR.REFT.REFE.  I PRT . ACCMT . ACCMB . MU . KPAR . KSTR . MVR . KK 
1 . JDEL , JDI S . MT . KT . MB . KB . RT . RB . BOB . BOT . FDS 
WRI T  E ( 50' 2 )  LDUMT . LDUMB . XDUMT , XDUMB , X 1 , X2 . LT . SF 
WRITE (50' 3)  PT.PB 

WRITE (50' 4)  KCONT , I MTS. LUNW , KDI S.DIST.KTYP. 1BTC, 1 5MP , XT  S , XBS , 

1 LT  S , SF  S , KPT , KPB . KIT, KIB. REST, RESB 
999  CONTINUE 

WR IT  E ( 5 1 ' KK ) I MT . XT , VLT . XB , BM . SM . FD . ( X 1 ( I ) . I  =  1 . 8 ) .  (  X2 ( I ) . 1  =  1 . 8 ) 

CALL  EXIT 

END 

//  FOR  PID 

SUBROUTINE  PID(KP.KI,RES. REF . XSP, USP , XTV . UCT , ACCUM) 

REAL  KP.KI 
UMI N= -O . 5 
UMAX =0. 5 
ER=(XTV-XSP)/REF 
RES  =  RES+  ER 
SP=KP*ER+KI *RES 
ACCUM  =  ACCUM  +  ER*  ER 
I F ( SP -UMAX )  1.1.2 
2  SP- UMAX 

GO  TO  3 

1  IF(SP-UMIN)  4,4,3 

SP=UMIN 
CONTINUE 

UCT  ~ ( 1 .0+SP)*USP 
RETURN 
END 
//  FOR  STR 

SUBROUTINE  ST R ( M . K . BO . R2 . P  .  X  .  LDUM . XDUM . XTN . USP . LTN . UCT . ACCUM . OUCT ) 
REAL  LTN, MU. LDUM( lO) 

DIMENSION  P(8.8).PHI(6.8),Y(8).X(8).PY(8), AKY (8,8), PPH( 8.8), 

1  E  RK ( 8 ) , XDUM( 8  ),AK(8) 

READ  (50'  1 )  KPT R. REFT .REFB,  I PRT, ACCMT. ACCMB. MU. KPAR.  KSTR, MVR. KK 
1 , JDEL 
DO  33  1 = 1 . 8 
DO  33  J= 1 , B 


' 


. 


. 


421 

r,H)(i.j)r0.o 

422 

33 

PHI(I.I)-i.0 

423 

UMAX® 1 . 5*USP 

424 

UMJfJ=0.  S^USP 

425 

L  +  K 

426 

n=l-*m 

427 

K2=K+2 

428 

K4=K+4 

429 

DO  204  1=1. K2 

430 

KIN=K2-I+2 

431 

204 

XDUM(KIN)=XDUM(KIN- 1  ) 

432 

DO  205  1=1 , K4 

433 

K I N=K4- I +2 

434 

205 

LDUM (KIN) =LDUM( KIN- 1 ) 

435 

DO  206  I = 1 , M 

436 

I K  1  =  I +K+ 1 

437 

206 

Y( I ) =XDUM( I K 1 ) 

438 

DO  207  I = 1 , L 

439 

IM=1*M 

440 

I K2= I +K+ 1 

441 

207 

Y ( IM) =LDUM( I K2 ) 

442 

XDUM ( 1 )=-XTN 

443 

LDUM ( 1 )=LTN*BO 

444 

IF(MVR-I)  50.51.50 

445 

50 

1 F ( KK-KSTR )  52.51.51 

446 

52 

CONTINUE 

447 

C  COMPUTE  THE  GAIN  VALUE 

448 

C 

449 

DO  22  1=1. N 

450 

PY ( 1  )  =0 . 0 

451 

DO  22  U= 1 , N 

452 

22 

PY(I )=PY(I)+P(I,J)*Y(J 

453 

YPY  =0 . 0 

454 

DO  23  1  =  1  ,N 

455 

AK ( I ) =FHI (1.1  ) *  PY ( I  ) 

456 

23 

YPY  =  YPY  + Y ( I )*PY( I  ) 

457 

YPY  =YPY+R2 

458 

DO  25  1=1. N 

459 

25 

AK( I ) =AK( I )/YPY 

4  60 

C  COMPUTE  COVARIANCE  MATR 

461 

C 

462 

DO  26  1= 1 ,N 

463 

DO  26  U=  1  ,  N 

464 

26 

AKY ( I , J )  =  PHI ( I . J)-AK( I 

465 

DO  27  I = 1 . N 

466 

DO  27  J=1 .N 

467 

PPH( 1 . J)=0.0 

468 

DO  27  I J= 1 , N 

469 

27 

PPH( I . J)=PPH(I , J)+P(I . 

4  70 

DO  24  1=1 ,N 

471 

DO  24  J= 1 , N 

472 

P(1 ,j)=o.o 

473 

DO  28  I U= 1 . N 

474 

28 

P( I . J)=P( I . U )  +  AKY ( I . I J 

475 

24 

P(I.d)=P(I.J)/MU 

476 

C  UPDATE  PARAMETER  VECTOR 

477 

C 

478 

YX=0 . O 

479 

DO  29  1=1. N 

480 

X(1 )=PH1 ( I . I )+X( I  ) 

9 


' 

' 
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<81 

29 

YX  =  YX-*  Y (  I  )  *X(  I  ) 

482 

C 

COMPUTE  ERROR 

483 

c 

Y ( 4  )  IS  EOUAL  TO  B0*U(T-2) 

484 

c 

485 

ERR=XTN-  LDUM(  K-*  1  )-YX 

48G 

DO  201  I-1.N 

487 

E  RK ( I )  =  AK ( 1 ) -ERR 

488 

201 

X (  I  )  =X (  I  )-*LRK(  I  ) 

489 

c 

COMPUTE  ACCUMALATED  ERROR 

4  90 

c 

491 

ACCUM  =  ACCUM+XDUM( 1 ) *  *2 

492 

c 

COMPUTE  CONTROL  VECTOR 

493 

c 

494 

51 

CONTINUE 

495 

YX=0.0 

496 

DO  202  I = 1 . M 

497 

202 

YX=YX*XDUM( I )*X( I ) 

498 

YY=0 . 0 

499 

DO  203  I  =  1  ,  L 

500 

IM=I+M 

501 

203  YY  =  YY+  LDUM( 1)*X(IM) 

502 

LTN=( -YY-YX )/B0 

503 

IF (KK-KPAR )  410.53,53 

504 

53 

CONTINUE 

505 

I F ( JDEL -  1 ) 99 , 99 , 98 

506 

98  UCT=OUCT+LTN*USP 

507 

GO  TO  97 

508 

99  UCT  = ( 1 . +  LTN ) *USP 

509 

C 

CONSTRAINTS  ON  THE  CONTROLLER 

510 

C 

511 

97 

I F ( UCT  -UMIN)44 ,44 ,45 

512 

45 

IF (UCT  -UMAX ) 1 IO, 47 , 47 

513 

47 

UCT  =UMAX 

514 

GO  TO  1 10 

515 

44 

UCT  =UM IN 

516 

GO  TO  1  10 

517 

4  10 

UCT  =  USP 

518 

1  io 

RETURN 

519 

END 

520 

CALL  GAUSS(IX.SA.AM.V) 

521 

AA=0. O 

522 

DO  7  1=1.12 

523 

I Y  =  I X*  899 

524 

IF  (IY)  5.6.6 

525 

5  IY=IY+32767+1 

526 

6  YFL=I Y 

527 

YFL=YFL/32767.0 

528 

I  X  =  I  Y 

529 

7  AA=AA+YFL 

530 

V=(AA-6.0)*SA+AM 

531 

RETURN 

532 

END 

533 

•DELET  A  RKT02  ***** 

534 

•STORECI I  A  1  RK702  RKT02 

535 

•FI LES( 50, DI STO.O) 

536 

•FILES(SI.DISTI.O) 

537 

•CCEND 

538 

//  END 

539 

END  OF  FILE 


0802 


020109 


. 

. 


. 


144 

APPENDIX  D:  TOP  PRODUCT  CONTROLLER  TUNING  BY  THE 
ZIEGLER-NICHOLS  CONTINUOUS  CYCLING  METHOD 


The  distillation  column  was  operated  at  the  following 
conditions: 


Feed  Flow 
Reflux  Flow 
Steam  Flow 
Top  Composition 
Bottom  Composition 
Feed  Composition 


18.00 

g/s 

13.70 

g/s 

14.70 

g/s 

96.00 

% 

0.00 

% 

50.00 

% 

Various  proportional  gains  were  tested,  under  closed 
loop  control,  for  a  1%  step  increase  in  set  point  to  obtain 
a  sustained  system  oscillation.  The  results  are  summarized 
in  Table  D 1 . 


Table  D1 

Results  of 

Top  Product  Composition 

Ziegler-Nichols  Tuning  Runs 

Run 

Propor  t iona 1 

Gain  Svstem  Performance 

D 1 

-18.0 

offset 

D2 

-12.0 

offset 

D3 

-12.0 

offset 

D4 

-20.0 

offset 

D5 

-30.0 

offset 

D6 

-50.0 

offset 

D7 

-75.0 

damped  oscillation 

D8 

-100.0 

sustained  oscillation 

The  control  and  manipulated  variable  responses  for  Run 
D8  are  plotted  in  Figures  D1  and  D2.  On  the  basis  of  the 
Ziegler-Nichol s  recommended  settings,  the  controller 
constants  for  proportional  plus  integral  control  are: 


Kp  =  0.45KPmaj:  =  -45.0 
TR  =  P*/1.2 


6.2  minutes 


b 


MEA5  (  PERCENT 
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TIME  CSEC] 


Figure  DI  Top  Product  Composition  Response  To  A  +1% 

Step  Change  In  Set  Point  Using  A  Proportional 
Controller  With  Kp  =  -100.0 


' 


REFLUX  (G/S) 


146 


Figure  D2  Reflux  Flow  Rate  Response  For  Proportional 

Controller  With  kp  =  -100.0  Of  Top  Product 
Composition  Control  For  A  +1%  Step  Change 
In  Set  Point 
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Appendix  E: 


Simulation  Program  Listing 


I 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 
1  1 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 
4  1 

42 

43 

44 

45 

46 

47 

48 
4  9 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 


4 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  •  c 

C  UNIVERSITY  OF  ALBERTA.  CHEMICAL  ENGINEERING  DEPARTMENT  C 

C  8  TRAYS  DISTILLATION  COLUMN  CONTROL  SIMULATION  C 

c  BY  RON  TONG  C 

C  c 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

r 


c 

c 

C 

C 

c 

c 

C 

c 

C 

c 

c 

C 

C 

c- 

C 

C 

c 

c- 

c 

c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  PROGRAM  USES  AN  EULER  METHOD  OF  INTEGRATION 
TO  PREDICT  THE  TRANSIENT  RESPONSE  OF  A  10  STAGE 
DISTILLATION  COLOUMN.  THE  ASSUMPTIONS  ARE 

1.  CONSTANT  WEIGHT  HOLDUP 

2.  UNIFORM  LIQUID  COMPOSITION-PERFECT  MIXING 

3.  NEGLIGIBLE  VAPOUR  HOLDUP 

THE  LIQUID  AND  VAPOUR  FLOWS  ARE  CALCULATED  IN  SUBROUTINE 
VF  LOW ,  USING  THE  TRAY  MASS  BALANCES  AND  ENTHALPY  BALANCES. 
THE  VAPOUR  COMPOSITIONS  ARE  CALCULATED  FRO  THE  LIQUID 
COMPOSITIONS  IN  SUBROUTINE  CLMP.  ENTHALPIES  FOR  LIQUID  AND 
VAPOUR  STREAMS  ARE  CALCULATED  IN  SUBROUTINE  ENTHALPY.  THE 
MAINLINE  PROGRAM  PERFORMS  THE  ACTUAL  I  IT  EGRAT ION  AND 
THE  SUBROUTINES  IN  THE  CORRECT  ORDER. 

CONTROL  STRATEGIES  AVAILABLE 

1.  PROPORTIONAL  INTEGRAL 

2.  SELF  TUNING  REGULATOR 


VARIABLE 

A I  E 
A I  EB 
AHLO 
EHLO 
BOB 
BOT 
D 
DA 
DT 
DXB 
FL 
HF 
HSI 
IB 
IBD 

ICHNG 

IN 

IR 

IV 

ISPC 

1ST 

ITEST 

UCONT 


UM 

KCONT 

KB 


FUNCTION 

SUM  OF  ERROR 
SUM  OF  ERROR 
INTERCEPT  IN 
SLOPE  IN 


S'"'  E  TERMS  IN 
SL  vE  TERMS  IN 
ST  E  M  ENTHALDHY 
STEAM  ENTHALPHY 


TOP  COMPOSITION 
BOTTOM  COMPOSITION 
RELATION 
RELATION 


CONTROL 
CONTROL 
TIME  (G/S) 


SCALING  FACTOR  '  SED  IN  STR  ’  OTTOM 
SCALING  FACTOR  USED  IN  STR  7  OP 
DISTURBANCE  VECTOR  AT  3  DIFFERENT 
FIRST  HALF  OF  DATE  TO  BE  PRINTED 
INTEGRATION  TIME  STEP  (SECOND) 

BOTTOM  COMPOSITION  DELAYED  BY  GC  MEASUREMENT  TIME  LAG 
INITIAL  FEED  FLOW  (G/S) 

FEED  ENTHALPHY  (KU/KG) 

STEAM  ENTHALPHY  (KJ/KG) 

BOTTOM  GC  MEASUREMENT  TIME  LAG  LOOP  COUNTER 
BOTTOM  G  C  MEASUREMENT  TIME  DELAY 
IN  NUMBER  OF  SAMPLING  INTERVAIS  (S.I.) 

VECTOR  OF  TIME  AT  WHICH  DISTURBANCE  OCCURS  (S.I.) 

INDEX  TO  CHECK  PRINT  TIME  FOR  COLUMN  DYNAMICS 

REFLUX  FLOW  NOISE  GAUSSIAN  CONSTANT  USED 

FEED  FLOW  NOISE  GAUSSIAN  CONSTANT  USED 

TIME  AT  WHICH  SET  POINT  CHANGES  ARE  TO  BE  MADE  (S.I.) 

INDEX  TO  CHECK  THE  SAMPLING  TIME 

TYPE  OF  DISTURBANCE 

CONTROL  CONFIGURATION  KEY 

UCONT  = 1  TOP  CONTR  L  UC0NT=2  BOTTOM  CONTROL 
UCONT  =  3  TOP  AND  EUTTOM  CONTROL 


PRINT  CONTROL  KEY 
KCONT  *  1  OPEN  LOOP 
PROCESS  TIME  DELAY 


(S.I.) 

KCONT  r  2  CLOSE  LOOP 
FOR  BOTTOM  CONTROL  (S 


I  .  ) 


* 


' 


. 


. 

. 


■ 

■ 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

1O0 

101 

102 

103 

104 

105 

106 

107 

108 

109 

1  lO 

1 1  1 

112 

113 

1 14 

1  15 

116 

117 

118 

1  19 

120 


148 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


KT 

KPB 

KPT 

KPAR 

KPIT 

KP1B 

KSTR 
KST 
I  KT 

LB 

LT 

LF 

LDUMB 

LDUMT 

LTS 

MB 

MT 

MU 

NI 

NTT 

NUMRUN 

OLT 

OSF 

PB 

PT 

i.P 

RB 

RT 

REFB 

REFT 

RESETB 
RESETT 
SO 
SDR 
SF 
SFS 
TE 
TR 
TS 
UA 
Uf'IN 
UMAX 
XBS 
XBSP 
XFT 
XTS 
XT  SP 
XOUMB 


PROCESS  TIME  DELAY  FOR  TOP 
PROPORTIONAL  GAIN  CONSTANT 
PROPORTIONAL  GAIN  CONSTANT 

time  when  str  is  effected 

PRIOR  TO  THIS 
BEYOND  STR  IS 
PRIOR  TO  THIS 
BEYOND  STR  IS 
TIME  WHEN  STR 
PROGRAM  S 
I KT  = 1  FOR 


) 


PI 

PI 


CONTROL  (SI 
USED  FOR  BOTTOM 
1  SED  FOR  TOP 
IS. I . ) 

TIME .TOP  SYSTEM  CONTROL  IS  ON 
EFFECTED  (S.I.> 

TIME. BOTTOM  SYSTEM  CONTROL  IS 
EFFECTED  (S.I.) 

PARAMETERS  ESTIMATION  IS  STOPPED 
MPLING  TH  E  COUNTER 
INITIATION,  I KT  =  3  AFTER, 


CONTROL 

CONTROL 


PI 


ON  PI 


(S.I.) 


FOR  USE  IN  SUE  OUT  I NE  COMP  ENTH  VFLOW 
NUMBER  OF  PARAMETERS  IN  MANIPULATED  VARIABLE 
FOR  BOTTOM  STR  (BETA) 

NUMBER  OF  PARAMETERS  IN  MANIPULATED  VARIABLE 
FOR  TOP  STR  (BETA) 

FEED  STAGE 

DUMMY  VECTOR  TO  S  ' ORE  PREVIOUS  DIFFERENCE  IN 
MANIPULATED  VARIABLE  FOR  BOTTOM  STR 
DUMMY  VECTOR  TO  STORE  PREVIOUS  DIFFERENCE  IN 
MANIPULATED  VARIABLE  FOR  TOP  STR 
INITIAL  REFLUX  FLOW  RATE  (G/S) 

NUMBER  OF  PARAMETERS  IN  OUTPUT  VARIABLE 
FOR  BOTTOM  STR  (ALPHA) 

NUMBER  C"  PARAMETER’S  IN  OUTPUT  VARIABLE 
FOR  TOP  STR  (ALPHA) 

EXPONENTIAL  FORGETTING  F ACT L A  FOR  STR  ALGORITHM 
SIMULATION  RUN  TIME  (MINUTE) 

NUMBER  OF  TRA-  IN  THE  COLUMN  +  REBOILER 
SIMULATION  RUN  NUMBER 

TOP  MANIPULATED  VARIABLE  REFLUX:  0LT=U(T-2)  IN  STR 
BOTTOM  MANIPULATED  VARIABLE  STEAM:  0SF=U(T-2)  IN  STR 


ALGORITHM 

ALGORITHM 

ALGORITHM 

ALGORITHM 


BOTTOM 

TOP 

BOTTOM 

TOP 


CONTROL 

CONTROL 

CONTROL 

CONTROL 


COVAR  ANCE  MA  RIX  FOR  STR 
COVAR’ANCE  MAiRIX  FOR  STR 
HEAT  .OAD 

NOISE  COVARIANCE  FOR  STR 
NOISE  COVARIANCE  FOR  STR 
REFERENCE  BOTTOM  COMPOSITION  USED  FOR  CALCU!  ATING 
NORMALIZED  BOTTOM  '"OMPOS I T I  ON  DEVIATION  ERROR 
REFERENCE  TOP  .LOMI  3ITI0N  USED  FOT  CALCULATING 
NORMALIZED  TOP  COMPOSITION  DEVIATION  ERROR 
RESET  TIME  FOR  BOTTOM  PI  CON  ROL  (MIN) 

RESET  TIME  FOR  TOP  PI  CONTROL  (MIN) 

STANDARD  DEVIATION  OF  FEED  FLOW  NOISE  (G/S) 

STANDARD  DEVIATION  OF  REFLUX  FLOW  NOISE  (G/S) 

INITIAL  STEAM  FLOW  (G/S) 

INITIAL  STEADY  STATE  STEAM  FLOW  (G/S) 

SECCUD  HAL-  OF  DATE  TO  3E  FRINTED 
REBOLER  TEM  ERATURE  (C) 

ST E AM  CHEST  TEMPERATURE  (C) 

HEAT  .RANSFER  COEFFICIENT  <  i/S) 

Ml  JIMJM  CONTROL  £R  CONSTRAINT 
MA \ I  MUM  CONTROLLER  CONSTRAINT 

INITIAL  STEADY  STATE  BOTTOM  COMPOSITION  ( WT  FN  CH30H ) 
NEW  DESIRED  STEADY  STATE  BOTTOM  COMPOSITION  (WT  FN  CH 
FEED  COMPOSITION  ( WT  FN  CH30H ) 

INITIAL  ST  ADY  STATE  TOP  COMPOSITION  ( WT  FN 
NEW  DESIRED  STEADY  STATE  TOP  COMPOSITION  (WT 


DH ) 


CM30H ) 

FN  CH30H) 


DUMMY  VE  CT 
IN  BOTTC  •: 


,R  TO 
STR 


STORE  PREVIOUS  OUTPUT  VARIABLE 


[*"  X 

.  ' 


• 

121 

122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 

161 

162 

163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 


149 


C 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


XDUMT  DUMMY  VECTOR  TO  STORE  PREVIOUS  OUTPUT  VARIABLE 

IN  TOP  STR 


REAL  LT ( 20 ) . LT S . KPT . KPB , MU . LDUMT ( 10).LDUMB(  10) , L  T  N , XDUMB (  10) 
DIMENSION  ISTRN(3),PT(7,7),PB(7.7) ,X1(7) . XB ( 7 ) . XDUMT ( 10) 
DIMENSION  D  (  4  )  ,  ICHNG(4  )  ,  RE  (  3  )  .  DF  (  4  ) 

COMMON  XT ( 20 ) , YT ( 20 )  , LT . VT ( 20) , HLT ( 20 ) . MVT ( 20 ) . HL ( 20) , E ( 20 ) , XF 
1(4) . KK , T  PROD , T , MTT , I KT , I 

COMMON  WT( 20) . SF .HSI , AHLO.BHLO.UA , T R , F ( 4 ) , QLP ( 20 ) . DT , LF . HF . OP . TS 
COMMON  KPAR.KSTR . MVR . UMI N , UMAX , KP I T ,KP1B 
DATA  ISRN/'RE  ' , ' FEED' , ' FC  '/ 

DATA  RE  ST, RE  SB, XI . XB/O . 0 . 0 . 0 , 7  *  0 . 0 . 7 *0 . 0/ 

DATA  PT , PB , LDUMT , LDUMB , XDUMT , XDUMB/ 4 9  *0.0,49*0.0,  10*0.  ,10*0.0, 

1 10*0. , 10*0.0/ 

DATA  KTOP. KBOT/1 . 1/ 


901  FORMAT (F5.2.2T5) 
900  FORMAT ( 10FS . 2 ) 

902  FORMAT ( 1015) 

903  FORMAT ( 1 JF8 . 4 ) 

904  FORMAT ( A4 ) 

905  FORMAT (8F 10. 5) 

906  FORMAT ( F8 . 2 . 8FS . 4 ) 
909  FOR  AT( 10F8.5) 


THE  FOLLOWIN'  'EAD  AND  WRITE  ST ' tMENTS  CCRRES  OND  T0  THE  DATA 
INPUT  REOU I R - ^  FOR  THE  SIMULAT1  N  MODEL 


INITI  L  VALUES 

RE  0( " . 999)NUMRUN 
399  FO  (15) 

WRI  6,998)  NUMRUN 

998  FOP  (»****»*************  SI MUL ATION  RUN  NUMB E R = ' ,15, 

\  *  *  *  *  •  •  ) 

REA. 15.997)  DA.TE 
997  FOR  A 7  (  2 A 4  ) 

WRITE (6. 996)  DA.TE 

996  FOR MAT (/,'*'. 3  2  X .  'DATE  PERFORMED  BY  R  TO  JG',.X,  '*',//. 
1 2A4 ,  10X.  ' • '  ) 


/  **««**+**** 


*  '  .  45X , 


r  PUT  DATA  L  LD  FOR  NL  COLUMN  MODEL 

RE AD( 5. 904)1  TEST 
READ ( 5 , 902 ) I CHNG 

WR I T  E ( 6 . 1  1 2  1  )  .«.,,...***•*...*«** 

1121  FORMAT  ('..•**♦»♦*•**»****************'***********'* 

. . . 

995  FORMAT (//////. 2X.  *  INPUT  DATA  USED  FOR  COLUMN  MODEL'/) 


. 


' 
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181 

182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 
211 
212 

213 

214 
2  15 
216 
217 
2  18 

219 

220 
221 
222 
223 
T24 
..'25 
226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 


C 

C 

WR I T  E ( 6 , 9  86 ) 

WRITE (6, 994) 

994  FORMAT (2X.  'IT  ESI') 

WRITE (6. 904) I  TEST 
WRITE  (6.986) 

V.  RI  TE  (  6  , 993  ) 

993  FORMAT ( 2X . ' ICHNG' ) 

WRITE (6 ,i 02) ICHNG 
WR  lTE(6.£ri6) 

RE  AD ( 5 , 903 )D 
WR 1 T  E ( 6 , 992 ) 

992  FORMAT ( 2X , 'D'  ) 

WR I T  E ( 6 . 903 )D 
WR ITE(6,986) 

READ(5.902)LF. IR 
WR I TE ( 6 , 99 1 ) 

991  F0RMAT(2X.  'LF' ,2X. 'IR') 

WRITE(6.902)LF.IR 
WRITE (6. 986) 

RE AD ( 5 , 901 ) DT , N I .NTT 
WR I T  E ( 6 , 990 ) 

990  FORMAT (2X.  '  DT  '  .  4X .  ' NI ' . 4X  .  'NTT'  ) 

WR I T  E ( 6 , 90 1 ) DT , N I .NTT 
WRITE(6 .986) 

M T  T  =  NTT+ 1 

READ  (5.900)  LT ( MTT ) . SD . SDR 
WRITE ( 6 , 989 ) 

989  FORMA  (2X,'LT(MTT)  S.D.  S.D.R  ') 

WRITE  . 6,'  00)LT(MTT),SD.SDR 
WRITE l 6 . 986 ) 

R  E  A  0 ( 5 . 900 ! F  L , XF  T , HF 
WRITE (6, 983) 

988  ORMAT ( 2 X ,  FL  XFT  HF .  '  ) 

WR  I T  E ( 6 , 9C 0)FL,XFT,HF 
WR ITE(6,9c  3) 

READ(5.900  )SF . HSI , BHLO , AHLO , UA , TR 
WRITE (6 . 987  ) 

987  FORMAT (2X  ,  '  SF  '  .  5X  ,  '  HS  I  ' . 8X ,  'AHLO' .4X.  'BHLO' . 2X . ' UA '  , 7X ,  ' TR ' ) 

WR  I  T  E ( 6 , 900 ) SF . HS I . AHLO . BHLO . UA . TR 
WR ITE(6,986) 

C 

C  IN.  UT  DA 'A  USED  FOR  STR  AND  PI  CONTRC _  ALGORITHM 

c 

READ( 5 , 903 )  LTS.XTS.XTSP.SFS,  BS.XBSP 
WRITE ( 6 , 985 ) 

985  FORMAT(////. 2X. 'CONTROL  ALGORITHM  INPUT  DATA') 

WRITE (6. 986) 

986  FORMAT  (  ' - - - - - - - - 

1 - '  ) 

WR I T  E ( 6 .984) 

984  FORMAT!  X  .  ' LTS '  . 5X . ' XTS ' . 4X .  ' XTSP '  . 5X .  ' SFS '  . 5X . ' XBS ' . 5X .  'XBSP'  ) 
WRITE  (6.903)  LTS . XT S . XT SP . SF S . XBS . XBSP 
WRITE (6. 986) 

READ  (5,903)  UM1N.UMAX 
WRI  TE ( 6 . 982  ) 

982  FOR  .T(2X,'MIN.  MAX.  CONTROLLER  CONTRAINT') 

WR 1  ; ( 6 , 903 )  UMIN.UMAX 

UMIiJ*  1  .  0-*  UM1  N 


X 


**..T 


■ 


. 


•• 


241 

242 

243 

244 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 

261 

262 

263 

264 

265 

266 

267 

268 

269 

270 

271 

272 

273 

274 

275 

276 

277 

278 

279 

280 

281 

232 

233 

234 

285 

286 

287 

288 

289 

290 

291 

292 

293 

294 

295 

296 

297 

298 

299 

300 
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UMAX* 1 .0+ UMAX 
W«!TE(6.986) 

READ  (5,902)  JM , <JC0NT .KCONT , KPAR ,KSTR ,MVR  KPIT  KPIB 
WHITE (6. 981  ) 

98  1  FORMAT  (  2X  ,  '  JU  UCONT  KCONT  KPAR  KSTR  MVR  KPIT  KPIB') 

WRITE  (6,902)  JM , UCONT , KCONT , KPAR . KST R . MVR , KP 1 T  KPIB 
WRITE(6.986) 

READ  (5.900)  KPT . KPB . RE SE TT . RE SE  1  B 
WRI T  E ( 6 , 980) 

980  FORMAT (IX, 'KPT  K  B  RESETT  RESETB') 

WRITE  (  6 , 900)  KPT . KPB . RLSETT , RESETB 
WRITE (6, 986) 

RE-'  0(5 , 902  )  MT,KT  .MB.KB.  IBD.  IV,  LTP,  LB 
WRITE (6. 979) 

979  FORMAT  ( 2X ,'  MT  KT  MB  KB  IBD  IV  LTP  LB') 

WRITE (6, 302)  MT.KT.MB.KB. IBD. IV. LTP. LB 
WR ITE(6,986) 

READ( 5 , 905 )  MU . BOT , BOE . RT . RB 
WRITE(6,978) 

978  F0Rf1T(2X,'  MU  BOT  BOE  RT  RB') 

WRIT £(6.905)  MU . BOT , BOB , RT  .  RB 
WRITE(G .986) 

READ  (5,900)  (PT(I,I),I  =  1,7),(PB(I  . I), 1  =  1. 7) 

WRIT E( 6, 977) 

977  FOR  AT ( 2X .  ' STR  XD  S  XB  INITIAL  CC7ARIANCE  MATRI >  USED') 

WRI  .  E(6 .900) (PT( I . I ) . 1= 1 ,7) . (PB( I , I ) , 1  =  1 ,7) 

WRITE (6. 986) 

RE AD (5. 909 )  (X1(I),I=1.7).(XB(I).I=1.7) 

WRI  :  ".(6,976) 

976  FOR  AT ( 2  X ,  'STR  XD  &  XB  INITIAL  P.  'AME  ERS  USED') 

V  -?]  T  E  (  6 , 909  )  (  X  1  (  I  ),I  =  1,7),(XB(I  )  1  =  1.7) 

WRITE (6. 986) 

C - 

c 

c 

c 

c 

c - 

c 

c 

C  USING  SUBROUTINE  INII  :  READS  IN  II.  TIAL  D I  ST  I  L  AT  I  ON  COLUMN 

C  VARIABLES-XT. YT.LT.VT.WT.OLT.E.HLT.HVT  FOR  EACH  TRAY  AVAILABLE 

C  FROM  DATA  FILE  LINE  22-35,  SUBTOUTINE  PRINT  :P  INTS  THEM 

C  USING  LUN=6 

F (  1)  =  FL 
X  F ( 1  )  =  XFT 
E  (  10)=0.0 
<  ALL  INII 
'  *  1 
1=0.0 
OP =327 4 7 . 

T  S= 107.7 
CALL  PRINT 

CONVERTING  FROM  NUMBER  OF  SAMPLING  INTERVALS  I  >0  NUMBER  OF 
INTEGRATION  TIME  STEPS  1  SAMPLING  INI ERVAL=64  SECONDS 
JM* JM" 64/DT 
KPIT=KPI T/DT-64 
‘vPIB  =  KPIB/DT*64 
..STR  =  KSTR/DT  *64 


' 


- 


■ 

- 


301 

302 

303 

304 

305 

306 

307 

308 

303 

3  10 

31  1 

312 

313 

314 

315 

316 

317 

318 

319 

320 

321 

322 

323 

324 

325 

326 

327 

328 

329 

330 

331 

332 

333 

334 

335 

336 

337 

338 

339 

340 

34  1 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 

353 

354 

355 

356 

357 

358 

359 

360 


ooo  oooo  ooo  ooo 


152 


KPAR=KPAR/DT»64 
NI =NI/DT*60 
KST  =64 /DT 
I SPC=20/DT»64 
I B  =  0 
I  ST  =0 
I N  =  0 

RE(  1  ) =LT ( MTT  ) 

I KT  =  1 


DETERMINING  THE  TYPE  OF  DISTURBANCE 
DO  5  1=1,3 

I F ( I STRN( I )-ITEST)5,7,5 

5  CONTINUE 

7  GO  TO  (6,8,9). I 

SETTING  UP  DISTURBANCE  VECTOR  ACCORDING  TO  TIME  VECTOR 

6  DO  61  1=1,3 

1 CHNG ( I ) = I CHNG ( I )/DT*64 
XF ( I )  —  XF (  1  ) 

DF ( I  )  =  F(  1  ) 

F(I)=F(1) 

61  RE ( I ) =D ( I  ) 

GO  TO  1 1 

8  DO  8 ’  1=1,3 

I CHN3 ( I ) = I CHNG( I )/DT*64 
XF( I  )  =  XF (  1  ) 

: E ( I  )  =  RE(  1  ) 

DF ( I ) =D( I  ) 

81  F ( I  ) =D ( I  ) 

GO  TO  11 

9  DO  91  1=1,3 

I CH'.'G (  I  )  =  IC  \'G(  I  )/D  *64 
RE ( I  )  =RE (  1  ) 

DF ( I  )  =  F  (  1  ) 

F(I)=F(1) 

91  XF( I  )=D( I  ) 

11  CONTINUE 

AT  THE  START  OF  RUN.  SUM  OF  ERROR  SQUARES  ( A I E  A I EB )  ARE  ^ET  TO  ZERO 
REFERENCE  TOP  AND  E  -TTOM  COMPOSITIONS  ARE  ALSO  DESIGNATED 

1  =  1 

A I E  0.0 
A I  E1S=0 . 0 
DXB=XT ( 1 ) 

XT  1  =T  .3 
REFT  =  . 96 
REFB  =  0 . 06 


Si;  JLATION  RUN  OUTPUT  LABEL 
WR I T  E ( 7 .  1 1 1 9 ) 

1119  F0RMAT(////.2X,  '  **********  SIMULATION  RESULTS  **********'.//.'  TI 
1ME  BOTTOM  STEAM  TOP  REFLUX  FEED  BOTTOM  ACCUM. ERROR 

2 SQUARE '  ,  /  ,  '  (MIN)  MEAS  (G/S)  (G/S)  ACTUA 

1 L  T  3'  ) 

WRITE (7 .906) r .XT 1.SF.XT(10).LT(10).F(I),XT(1).AIE.AIEB 


. 
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361 

362 

363 

364 

365 

366 

367 

368 

369 

370 

371 

372 

373 

374 

375 

376 

377 

378 

379 

380 

381 

382 

383 

384 

385 

386 

387 

388 

389 

390 

391 

392 

393 

394 

395 

396 

397 

398 

399 

400 

401 

402 
40  > 
40 

405 

406 

407 

403 
4o: 
4  10 
411 
4  12 
4  13 
414 
4  15 
4  16 
417 
4  18 
4  19 
420 


WB  J  HE ( 8 . 1 122) 

1122  FORMAT*//////'*  *******  PARAMETERS  VALUES  ♦•••••***•'/) 

'.'RITE  (8. 1  123) 

1123  DRMAT ( 8X , ' S  T  R  CONTROL') 

WRI T E ( 8 . 1124)  MT,  KT .  LTP 

1124  FORMAT ( 8X ,  ' M* " , I  2 , 2X ,  ' K* ' , 1 2 , 2X , ' L= '  12/) 

WRJTE(8. 1 125) 

1125  FORMAT ( 3X  .  ' TIME '  . 2X .  ' ALPHA 2 '  , 2X ,  ' BE  T  A2 '  , 4X ,  'ALPHA  1 ' ,2X. 'BETA  1 
12X..  'BETA3'  ) 

C 

C 

c 

c 

c 

c 

C  COLUMN  SIMULA) ION  ITERATION  LOOP  STARTS 

C 

c 

c 

DO  lOO  KK  = 1 ,NI 
T-KK*  DT / 60 . 

C 

C  0  CHECK  TIME  FOR  SET  POINT  CHANGES 

K(KK.GE - ISPC)  XTS  =  XTSP 
I F ( KK . GE . I SPC )  XBS=XBSP 
AIE  =  AIE  +  (XTS-XT( 10) )*  f  *  . S-XT(  10) ) 

AIEB~AIEB+(XBS~XT ( 1 ) )  (XBS-XT(I)) 

c 

C  TO  C  iECK  TIME  FOR  LOAD  DISTURBANCE 

DO  23  JJ=  1 , 3 
II  "<K  ICHNG(JJ)  >23,22.23 

22  I--JJ 

23  CONTINUE 
C 

C  03T;  N  CORRECT  FEED  AND  TOP  PR-DUCT  FLOW  RATES. THEN  SIMULATE 

FFF=D<  (  I  ) 

25  i PROD= VT (  r)-LT(MTT) 

C 

c 

C  SUB NO  /TINE  COMP:  SING  TABLE  LOOK  UP. WITH  LIQUID  COMPOSITION 

C  (XT)  AVAILABLE.  FiND  F  JILIBRIUM  VAPOR  COMPOS ITIC.J(YST); 

C  AND  WITH  EFFICIENCIES  SS I GNED . AlTUAL  '  \POR  COMPOS  I T 1  ON ( Y T ) 

C  FOR  EACH  TRAY  IS  OBTAj  3D 

CALL  COMP 

c 

c 

C  su:  OUT  1 NE  ENTH:  USING  TABLE  LOOK  yP  .LIQUID  ENTHALPHY ( HL  F ) 

C  AND  VAPOR  ENTHALPHY  ( T  ,'T  )  CAN  BE  OBTAINED  FROM  XT.YT 

10  CALL  ENTH 
C 
C 

C  SUBRC  TINE  VFL'W:  LOWS  LT  VT  RE  CALCULATED  USING  TRAY  MASS 

C  AND  E  THALPY  E  LANCES 

CALL  VF  LOW 
107  IKT  =3 
C 

c 

C  SUBROUTINE  INTGR  PERFO  \MS  ACTUAL  INTEGRATION 

c 


CALL  INTGR 
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421 

C 

422 

C 

CHECK  PRINT  TIME 

423 

JN= IN+ 1 

424 

I F ( IN- UM  )  20 , 52 , 20  , 

425 

52 

CALL  PRINT 

426 

1*4  =  0  • 

427 

20 

CONTINUE 

428 

c 

429 

c 

CHECK  SAMPLING  INTERVAL 

4  30 

104 

I ST  =  I ST+  1 

431 

IF  (1ST-KST)  100,883,100 

432 

888 

I  ST=0 

433 

c 

ROUTINE  GAUSS  IS  USED  TO  GENERA ! E  NORMALLY  DISTRIBUTED 

434 

c 

RANDOM  NUMBER  WITH  STANDARD  DEVIATION  0.05  G/S 

435 

c 

FOR  ADDING  NOISE  TO  THE  FEED  FLOW 

436 

CALL  GAUSS( IV.SD.O.O. VF ) 

437 

c 

REFLUX  FLOW  NOISE  CAN  BE  ADDED  FOR  OPEN  LOOP  IF  REQUIRED 

438 

c 

BY  DEFINING  SDR  TO  BE  NON-ZERO  QUANTITY 

439 

CALL  GAUSS( IR.SDR.O.O.VR) 

440 

F ( 1 )=FFF+VF 

441 

IF  (KCONT-1)  369,369,963 

442 

369 

LT  MTT )=RE( I  ) 

443 

LT ( MTT ) =LT ( MTT )+VR 

444 

c 

CHECK  FOP  OPEN  LOOP  OR  CLOSED  LOOP  SIMULATION 

445 

c 

446 

963 

GO  TO  (222,666),  KCONT 

447 

666 

CONTINUE 

448 

c 

CHECK  FOR  PI  OR  SIR  CONTROL - 

449 

c 

4  50 

I  F  (  KK  KP  IT)  53,53,57 

451 

c 

CHECK  FOR  TOP  CONT  OL .  BOTTOMS  CONTROL  SR  SIMULTANEOUS  CONT 

452 

c 

453 

53 

GO  TO  (33,57,33) . UCONT 

454 

33 

CALL  PID  KPT, RESETT, REST, XTS.LTS,XT( 10). LT (10), REFT) 

455 

OLT=LT ( 10) 

456 

57 

I F ( KK  KP IB)  56,56,777 

457 

56 

GO  TC  (40,44,44).  UCONT 

458 

44 

IB=1B+  1 

459 

c 

CHECK  BOTTOM  G.C.  MEASURE  ENT  TIME  DELAY 

460 

c 

461 

IF(IB-IBD)  40.41.40 

462 

4  1 

I B  =0 

46 

XT  1=DXB 

464 

DXB  =  XT (  1  ) 

465 

CALL  PI0(KPB.RESETB.RESB.XES.SFS.XT1 . SF , REFB ) 

466 

CSF  =  SF 

467 

40 

CONTINUE 

468 

777 

CONTINUE 

469 

IF  (KK  .  GE .  KP'R  .AMD.  KK  .LE.  KPIT)  GO  TO  ^2 

470 

IF(KK-KPIT)  63.63.62 

471 

62 

GO  TO  (55,63.55).  UCONT 

472 

c 

COMPUTE  THE  NORMALIZED  DEVIATION  VARIABLES 

473 

c 

474 

55 

I F ( KT OP . EQ .  1 )  OLT  =  L  T (  0) 

475 

KT  OP  =  2 

476 

XTN=(XT( 10) -XT S) /REFT 

477 

LTN= ( LT ( 10)-0LT )/LTS 

478 

OLT  T  =0L  T 

479 

OLT  =  LT (10) 

480 

CALL  s T R ( M T . KT . LTP . MU  BOT . RT . PT . X 1 . LDUMT . XDUMT .XTN.LTN.LTS. 

. 


- 
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481 

1LT( 10) , OLT ) 

482 

WRITE  (8.903)  T.X1(2).X1(4).X1(1).X1(3).X1(5).X1(6),X1(7) 

483 

63 

IF(KK-KPIB)  222.222.65 

484 

65 

I  F  (  KBOT  .  EQ.  1  )  OSF--SF 

485 

KBOT  =  2 

486 

GO  TO  (42.66.66).  JCONT 

487 

66 

I B= IB  + 1 

488 

IF(IB-IBD)  42.43.42 

489 

43 

I E  =0 

490 

XT  1 =DXB 

431 

DXB=XT ( 1 ) 

492 

C 

COMPUTE  NORMALIZED  DEVIATION  VARIABLES 

493 

C 

494 

XEN= ( XT 1-XBS)/REFB 

495 

SFN=(Sr-OSF )/SFS 

496 

OSFF  =0  jF 

497 

OSF  =SF 

498 

.<BN  =  -XBN 

499 

CALL  STR ( MB . KB . LB . MU , BOB . RB . PB . XB . LDUMB . XDUMB . XBN. SFN . SFS . SF . OSF  F ) 

500 

WRITE  (9.906 )T . XB ( 2  )  . XB ( 4  )  . XB (  1  )  .  XB ( 3 ) . XB ( 5 ) . XB ( 6 ) . XB ( 7 ) 

501 

42 

CONTINUE 

502 

222 

WRITE (7 ,906)  T . XT  1  . SF . XT ( 10 ) . LT ( 10 ) . F ( I ) . XT ( 1 ) . A I  E  .  A 1  EB 

503 

W  R I T  E (  1 1 .531 )XT( 10) ,LT(  10) 

504 

531 

FORMAT (2F8. 4) 

505 

lOO 

CONTINUE 

506 

WR I TE ( 7 , 905 )  T.AIE.AIEB 

507 

STOP 

508 

END 

509 

C 

510 

c 

511 

c 

SELF-TUNING  REGULATOR  CONTROL  ALGORITHM 

512 

c 

513 

SUBROUTINE  STR  (  M .  K  .  L  .  MU .  BO .  R2  .  P  .  X  .  LDUM  .  XDUM .  XTN' .  LTN.USP.UCT.OUCT  ) 

514 

REAL  LT  (  20 )  ,  f  "J  ,  LTS  .  LTN ,  LDUM(  1 0 )  ,  LTGS  .  XDUM(  10 )  ,  K(7) 

515 

DIMENSION  P(7,7).PHI(7,7),Y(7).X(7).PY(7).AK(7).AKY(7,7),FFH(7,7) 

516 

COMMON  XT (20) . YT(20) . LT , VT ( 20) , HLT ( 20) , HVT ( 20) . HL ( 20) . E  (  20) .XF 

517 

1(4) , KK , TP ROD , T . MTT , IKT . I I 

518 

COMMON  WT( 20) . SF .HSI . AHLO .EHLO.UA . ; R . F ( 4 ) . QLP ( 20 ) , DT . LF . HF , QP . TS 

519 

COMMON  KPAR . KST  R , MVR . UMIN . UMAX ,KP1T ,KPTB 

520 

DATA  LUNR.LUNV . ACCUM . PHI /5 . 6 . 0 . 0 , 49 *0 . 0/ 

521 

c 

•  •  , 

•4' 

=  NO.  OF  PARAMETERS  IN  Y.  K=  NO.  OF  TIME  DELAYS 

522 

c 

523 

ULCW=UMI N*USP 

524 

UHI =UMAX *  USP 

525 

N  =  L  +  M 

526 

K2=K+4 

527 

K4=K+6 

528 

0  208  1=1. N 

529 

208 

.  HI ( I . I  )=  1  .0 

530 

c 

GET  THE  DATA  FROM  PROCESS  MODEL 

531 

c 

532 

DO  204  I  =  1  , K2 

533 

KlN=K2-I+2 

534 

T4 

XDUM( KIN) =XDUM( KIN- 1 ) 

535 

DO  205  1=1, K4 

536 

K I  N  =  K4  -  I  2 

537 

05 

LDUM ( K I N )  =  LDUM( K I N- 1 ) 

538 

DO  206  I = 1 . M 

539 

IK1  =  I-*K+1 

540 

206 

Y ( I )=XDUM( I K 1 ) 

■ 
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54  1 

DO  207  I-1.L 

542 

I M  =  I  +  M 

543 

I  K2  =  1  -*  K-M 

544 

207 

Y( IM)=LDbM( IK2)  ' 

545 

XDUM (  1  )=-XTN 

546 

LDUM  (  1  )=LTN*B0 

547 

IF(MVR-I)  50.51.50 

548 

50 

I  F  (  KK-KST  R  )  52.51,51 

549 

52 

CONTINUE 

550 

C 

COMPUTE  THE  GAIN  VALUE 

551 

C 

552 

DO  22  1=1, N 

553 

PY( I )=0.0 

554 

DO  22  J= 1 , N 

555 

22 

PY(I)=PY(I)4P(i,j)*y(J) 

556 

Y  P  Y  =  0 . 0 

557 

DO  23  1=1, N 

558 

AK ( I ) =PHI ( I . 1 ) ' PY ( I ) 

559  . 

23 

YPY=YPY+Y ( I )*PY( I ) 

560 

YPY  = YPY+R2 

561 

DO  25  I = 1 . N 

562 

25 

AK( I  ) =  AK( I  ) / Y  P Y 

563 

C 

Compute  covariance  matrix 

564 

C 

565 

DO  26  1=1. N 

566 

DO  25  J= 1 , N 

567 

26 

AKY ( I  , U ) =PHI ( I  . J ) - AK  f I ) * 

568 

DO  27  I =  1  . N 

569 

DO  27  0=1. N 

570 

PPH( I . J)=0.0 

57  1 

DO  27  I 0= 1 , N 

572 

27  PPH(  I  ,  U  )  =  PPH(  I  .  ;)  +  p(i.i-- 

573 

DO  24  1 =  1  ,N 

574 

DO  24  0;  1  ,N 

575 

P( I . j)-o.o 

576 

DO  2  3  I J  = 1  .  •  J 

577 

28 

?(I  .  0)=P(  1  .  d)->  AKY<  I  .  10  * 

578 

24 

5 ( I , 0 ) = ( P ( I . 0 ) ) /MU 

579 

C 

UP i.  ATE  PARAMETER  VECTOR 

580 

C 

581 

YX=0 . 0 

582 

DO  29  I = 1 . N 

583 

X(I  ) =  PHI ( I  .  I)’X(I) 

584 

29 

YX  =  YX  +  Y ( I  )  <(I) 

585 

C 

COMPUTE  ERRC 

566 

C 

Y ( 4 )  IS  EQUAL  TO  B0*U(T-2) 

587 

C 

588 

err=x :n-ldu  (K+1 )-YX 

589 

DO  201  1= 1 .  J 

590 

ERr;(  I  )  =  A K (  I  )  ’ERR 

91 

201 

X( I  )=X( I  )  +  E RK ( I ) 

92 

C 

(  JMPUTE  CONTROL  VECTOR 

593 

c 

594 

51 

CONTINUE 

595 

YX=0 . 0 

596 

DO  202  I-1.M 

597 

202 

YX  =  YX  +  XC  *U)‘X(I) 

598 

Y  Y  •  T.O 

599 

DO  203  ‘  1 . L 

600 

1 M  =  I  +  M 

PPH( I J. J) 


' 


■ 


601 

602 

603 

604 

605 

606 

607 

608 

609 

610 

611 

612 

613 

614 

615 

616 

617 

618 

619 

620 

621 

622 

623 

624 

625 

626 

627 

628 

629 

630 

631 

632 

633 

634 

635 

636 

637 

638 

639 

640 

64  1 

642 

64  3 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

1 1 

12 

13 

14 

15 

16 

17 
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203  YY  =  YY  +  L  DUM (  1  +0  )  ‘  X  (  I  f.*. ) 

LTN=(-YY-YX)/B0 
IF(KK-KPIT)  46.46.53 
53  CONTINUE 

C  CONSTRAINTS  ON  THE  CONTROLLER 

C  * 

UCT  =OUCT  +  LTN*USP 
1F(UCT-Ul 0W)44 .44 .45 
45  IF(UCT-UHI >46,47.47 

47  UCT  =  UHI 

CO  TO  46 
44  'JCT  =  ULOW 

46  CONTINUE 

WR1 TE ( !  UNW . 37 )  ( Y ( I ) . X ( I  )  . P ( I . I ) . AK( I ) . I  =  1 . N ) 

37  FORMAT (4(7X,E10.3)) 

40  CONTINUE 

RETURN 
END 
C 
C 

C  PROPORTIONAL  INTEGRAL  CON  OL  ALGORITHM 

C 

c 

SUBROUTINE  PID( KP , RESET . RES. XSP ,USP , XT V. UCT . REF ) 

REAL  KP. LT( 20) . LTS . LTGS 

CCV.MON  XT  (20)  ,  YT  (  20 )  .LT.VT  (20)  ,HLT(20)  ,  HVT  (  20 )  .  HL  (  20)  .  E  (  20)  .  XF 
1(4)  .KK.TPR'  0  .  T  ,  !-1TT  ,  I  KT  ,  I 

COMMON  T  (  )  .  SF  ,HSI  .  AHLR  ,  BHLr' .  UA  ,  TR  ,  F  (  4  )  ,  QLP  (  20  )  .  DT  .  LF  ,  HF  .  OP  .  TS 

COMMON  K:  . KSTR . MVR , UMI N , UMAa 

UL0W=UM1'  USP 

UHI =  UMA  .  USP 

ER=(XTV-XSP)/REF 

RES=i  ;S+ER 

SP=Kr  ER+KP/RESET*RES 
UCT  =  (  1 .  +  SP)‘USP 
I F (UtT-UHI )  1.1.2 

2  UCT  =  UHI 
GO  TO  3 

1  IF(UCT-ULOW)  4.4.3 

4  UCT  =  UL  GW 

3  CONTINUE 
RETURN 
END 

■UC-OUTINE  GAUSS(IX.SA.AM.V) 

A  =  0 . 0 

DO  7  1=1.12 
I Y  = I X ’65539 
I  F ( I Y )  5.6.6 

5  IY=IY+2147483C i7+1 

6  YFI  Y 

YFl  :FL*0.465£j13E”9 
IX=  1Y 

7  AA=AA  YFL 

V=( AA-6 .0) ' SA+AM 

RETURN 

END 

SUBROUTINE  -RINT 

CCMM0NTXT?2O).YT(2O).LT.Vi  ( 20)  . HLT ( 20) . HVT ( 20) . HL ( .  0).E(20),XF 
1(4), KK . T PROD , T . MTT , IKT . I 


. 
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19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

4  1 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

7  1 

72 

73 

74 

75 

76 

77 
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COMMON  WT ( 20) . SF ,HS1 . AHLO . EHLO , UA . TR , F ( 4 ) . OLP( 20) . DT . LF . HF . OP . T S 
COMMON  KPAR  .  KST  R  .  MVR  .  UMIN .  UMAX  .  KP  I  T  .  KP  I  B 

101  FORMAT ( '  ' .20X. 'TIME ' .Ff  2. '(MIN)') 

102  FORMAT ( / . 2  X .  'STAGE'  .2X.  'LCOMP'  .2X.  ' VCOMP '  . 2X .  ' EF  FNCY ' .2X.'LFLOW',2 

IX. * VFLOW' , 2X . 'HOLDUP ' . 2X . ' LENTH ' , 2X . ' VENTH' )  , 

103  FORMAT ( / , 16.F7.3.F7.3.F8.3.F8.3.F7.3.F7.2.F9.2.F8.2) 

104  F  CKM.AT  (  /  .  5X  .  'FEED  RAT  E  '  .  2X  .  '  F  E  ED  COMP  '  .  2X  .  '  ST  E  AM  FL0W'.2X.'HT  COEF 
IF' „2X, 'HEAT  LOAD' ) 

105  FORMAT ( /. F 11.2.F11.4.F12.2.F12.3.F15.2) 

106  FORM -T ( / ,  10X .  ' RE BO  I LER  TEMP'.F9  3. 5X. 'STEAM  CHEST  TEMP'.F9.3) 

WRITE ( 6 ,  101  )  T 

VRITE(6. 102) 

DO  1  N=1  MTT 

1  WRITE (6.  03 )N . XT ( N) . YT ( N) , E ( N) , LT ( N) , VT ( N ) , WT (N) , HLT (N) . HVT ( N ) 
WRITE (6. 104 ) 

WRITE (6. 105)  F(I).XF(I),SF,UA,QP 

WRITE  (  6 , 1 06  )  T  R  .  T  S 

RETURN 

END 

SUBROUTINE  INTGR 
REAL  LT ( 20 ) 

DIMENSION  C  20) . FF (20) .XX( 20) . 0(20) 

COMMON  XT(20).YT( 20 ) . LT , VT ( 20)  , HLT(20),HVT( 20) , HL ( 20)  E ( 20 ) . Xr 
1 (4) ,KK, TPROD, T . MTT , I KT , I 

COMMON  WT ( 20 ) . SF . HS I . AHLO . BHLO . U A . T R . F ( 4 ) ,QLP(20),DT.LF.HF.QP.TS 
COMMON  KPAR , KSTR , MV R , UMI N . UMAX , KPIT.KPIB 
DO  25  N= 1 . MTT 

25  XX( N )  =  XT ( N ) 

NN=2 

26  DO  10  N*1,MTT 
K=N-  1 


J  =  N  +  1 

I F ( N- 1)1.4. 1 

1  1 F ( N-LF )2 , 5 . 2 

2  IF (N- MTT )3 . 6. 3 

3  D(N)  =  LT(J)*XT(J)+VT(K)*YT(K)-  N)  *  XT  (IJ  )  -  VT  (N )  *  YT  (  N) 

GO  TO  10 

4  D(N)=LT(J)*XT(J)-LT(N)’*XT(N)-VT(  Y  T  (  N  ) 

5  D(N)*LT?J)*XT(J)+VT(K)*YT(K)-LT(N)*XT(N)-VT(N)*YT(N)+F(I )*XF(I) 

GO  TO  lO 

6  D ( N ) =  VT (K)*(YT(K)-XT(N) ) 

10  CONTINUE 

I F ( NN - 3 ) 11.27,11 

11  DO  15  N=1 , MTT 
15  Q(N)=DT/W7(N) 

DO  22  N-  1 , MTT 
FF(N)=D(N) 

XT(N)=XT (N)+0(N)*D(N) 

22  CONTINUE 

1 (NN-2  )  27 , 28 . 27 


28  .'  =  3 

v  LL  COMP 
CALL  ENTH 
CALL  '.FLOW 
GO  TO  26 
27  CONTINUE 

DO  30  N= 1 .MTT 
30  XT(N)  -XX(N)+(Q(N)/2) 
RETURN 


(  D  ( N )  ♦  F  F  ( N )  ) 
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1  12 
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115 
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END 

SUBROUTINE  INI  I 
REAL  LT ( 20 ) 

COMMON  XT (20) , YT ( 20) . LT . VT ( .  0) . HLT ( 20) , HVT ( 20) . ML ( 20 ) . E ( 20 ) . XF 
1(4), KK .  I  PROD , T , MT  T , I KT . I 

COMMON  WT ( 20) , SF .HSI . AMLO.EMi  O.UA ,TR , F(4 ) ,OLP( 20) ,DT , LF ,HF ,6p. TS 
COMMON  KPAR.KSTR. MVR , UMIN , UMAX , KPI T , KPI8 
NTT  =  MTT -  1 

READ( 5 . 100) ( XT(N) ,  N  = 1 , MTT ) 

RE  AD  (  5 , 1 00  )(YT(N).N=1,’.TT) 

R  E  A  D  (  5 , 1 00  )(LT(N)  ,f.'=1  ,NTT  ) 

READ(5, 100)(VT(N),f.=  1  ,NTT) 

READ( 5 . 100) ( WT (N) , N= 1 , MTT ) 

READ( 5 . 100) ( CLT  IN). N= 1 , MTT ) 

RE AD( 5 , 100 ) ( r ( N ) . N- 1 . NTT ) 

READ  (  5 , 100)  ( HLT  (N)  .’.'=1  .MTT) 

READ( 5 . 101 ) ( HVT ( N ) ,N=1 , MTT ) 

100  FORMAT (  1 0F 8.4) 

101  FORMAT ( 10F  8 . 3 ) 

RETURN 
END 

SUBROUTI  f.'E  VFLOW 
REAL  LT(.O) 

DIMENSION  C( IE  ,  B ( 1 9 ) . F T ( 20 ) . SB ( 1 9 ) . SFT ( 1 9 ) .DHL (20) ,SW( 19) 
COMMON  XT (20) . V , 1 20) . LT. VT( 20) . HLT (20) . HVT ( 20 ) . HL ( 20 ) . E ( 20 ) . XF 
1  (  4 ) . KK . TPROD , T . MT T , I KT . I 

COMMON  WT( 20) . SF .HSI  . AHLO.BHLO.UA . TR . F  < 4 ) . OLP( 20) . DT , LF . HF . OP . TS 
COMMON  KPAR.KSTR. MVR. UMIN .UMAX.!  : I T . KP  I  B 
NTT  =  M7 T- 1 

-TEMPERATURE  OF  STEAM 
C* -HE  AT  ADDED 

I F ( I KT -2)60.51,51 
60  READ( 5 . 100) I QP 
OP=IQP 

: F( I  OP -0)50.5 •  , 50 
lOO  ORMAT (217 ) 

51  TS  =  ( SF  *  ( HSI  -  A.  -:L0  )  JA  -  TR  )  /  '  U  A  +  EHLO  *  SF  ) 

OR=UA* (TS-TR ) 

OP  =QR 
50  CONTINUE 
NL=NTT- 1 

C  C--  VAPOUR  AND  LIQ1  J  FLOURAETS 
DO  1  N= 1 , NL 
M  =  N+  1 

C ( N ) = ( HLT ( M ) -H\  (N)) 

1  B ( N ) = ( HLT ( N) -HLT ( M ) ) 

C ( NTT ) = - ( HVT ( NTT ) -HLT (NTT) ) 

c  c-**  der:vatives  of  hl 

DO  35  N= 1 .  TT 

35  DHL ( N ) = ( HL  N)-HL(N))/DT 
M = NTT  + 1 
'0  2  N=1 .NTT 

2  SW(N)=0.  ) 

DO  4  N -  .NL 
1F(N-  J.5.20 

20  IF(N-LF )21 .6.21 

21  1F(N-LF  (22.22.7 

22  FT(N)=W;  ( N ) -DHL (  • )  -  B ’  J ) * ( SW( N ) -F ( I ) -LT ( M ) )*OLP(N) 

GO  TO  4  ... 

5  FT  (  N  )  =\/T  (N)  *1  L('  -B(N)*(S\.'(N)-F(  I  )  -  L  T  ( M  )  )  •  0  L  P  ( N )  *  QP 
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60  TO  4 

S  FT (N) =WT(N)»DHL(N) -B(N)» ( SW ( N ) - L T ( M ) ) - F ( I ) * (HF-HLT (N) )+QLP(N) 

GO  TO  4 

T  FT(N)=WT(rj)*DHL(N)-B(N)*(SW(N)-LT(M))+QLP(N) 

*  CONTINUE 

FT  (NT  T ) =WT (NTT ) ‘DHL (NT T ) - LT ( M ) * ( HLT ( M ) -HLT ( NTT ) )+QLP(NTT  )' 

ADD^O.O 

DO  lO  N=  1  .NTT 

SFT(N)  =  ADD-»FT(N) 

Kt>  ADD=SFT(N) 

ADD  =  0  O 

DO  11  N=1,NL 

SB(N)=ADD+B(N) 

1*  ADD  =  SB ( N ) 

VT(  NTT  )  =  SFT  (  NTT  )  /  (C(NTT  )•*  SB(NL)  ) 

DO  12  N= 1 , NL 

*2  VT(N)=(SFT(N)-SE(N)*VT(NTT) )/C(N) 

DO  13  N= 1 , NTT 
K  =  NTT  + 1 -N 
U=K-1 
M  =  K-M 

IF  (  K-LF  )40 .  U.40 
-30  I  F  ( N-NTT  )  4  1  .  .5.4  1 
4f  LT(K)=VT ( J  )  +  LT (M)-VT(K) 

GO  TO  13 

14  LT(K)=VT( J)4LT(M)-VT(K)+F( I ) 

GO  TO  13 

15  LT (K ) =LT (M  - VT(K) 

'.3  CONTINUE 

RETURN 

END 

SUBROUTINE  i  NTH 
REAL  LT ( 20 ) 

DIMENSION  XST( 19) . YST( 19) ,HXS( 19) . HYS( 19) 

COMMON  XT (20) . YT(20) . LT,VT(20) , HLT ( 20 ) , HVT  20 ) , HL ( 20 ) . E ( "O) . XF 
1(4  )  .KK,  TPROD,  T  .1  »  T  .  I KT  ,  I 

COMMON  WT( 20) , SF ,HSI . AHLO .BHLO.UA . TR . F (  ’  ) . QLP( 20) . DT . LF .HF .  jP . TS 
COMMON  KP AR . KSTR . MV R . UMIN . UMAX , KP 1 T , KP I  3 

100  FORMAT ( 10F8 . 4 ) 

101  FORMAT (  10F8 .3) 

102  FORMAT ( 10F8 .2) 

I F ( I KT - 2  )  1 , 4 , 4 

1  READ ( 5 , 100)XST 
READ( 5 .  ■ 00 ) Y  ST 
READ ( 5 ,  ,01  )HXS 
READ( 5 , 102 )HYS 
DO  50  N= 1 . 19 

50  CONTI  fJUE 

4  DO  5  N  = 1 , MT T 

5  HL  (  N )  -  HLT  ( N  ) 

K=  1 

DO  10  N= 1 , MT  T 
DO  2C  J=K. 19 
1 F  f  .  T( J)-XT(N)  20.20.2 
20  COi.i  INUE 

2  K=J-1 

HLT(N)=((XT(N)-XST(K))*(HXS(L)-HXS(K))/(XST(L)-XST(K))4HXS(K)) 

10  CONTINUE 
K*1 
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198 

DO  30  N-  1  .(ITT 

199 

DO  40  J  =  K,  19 

200 

IF(YST( J)-YT(N) )40.40.3 

201 

40 

CONTINUE  '  » 

202 

3 

K-U-  1 

203 

L  =  il 

204 

HVT(N)=( (YT(N)-YST(K))*(HYS(L)-HYS(K))/(YST(L)-YST(K))+HYS(K) ) 

205 

30 

CONTINUE 

206 

RETURN 

207 

CND 

208 

SUBROUTINE  COMP 

209 

C 

C 

YST  IS  THE  VAPOUR  COMPOSITION  WHICH  WOULD 

210 

C 

C 

BE  IN  EQUILIBRIUM  WITH  THE  LIQUID  ON  THE 

211 

c 

c 

CORRESPONDING  TRAY.  THEY  ARE  CALCULATED 

212 

c 

c 

FROM  X-Y  DATA.  AND  ARE  EXPRESSED  AS  LINEARIZED 

213 

c 

c 

EQUATIONS.  PLATE  EFFICIENCIES  ARE  ALSO 

214 

c 

c 

CALCULATED  He RE . 

215 

REAL  LT ( 20 ) 

216 

DIMENSION  YST (20) ,YS(26) ,XST(26) 

217 

COMMON  XT (20) . YT(20) ,LT.VT(20) .  LT(20) ,HVT(20) . HL ( 20) . E ( 20) . XF 

2  18 

1(4)  ,  KK  ,  TPROD  ,  T  ,  MTT..  IKT  .  I 

219 

COMMON  WT  (  20 ) , SF , K3 I ,  AHLO.BHLO  . ’JA  .  TR  ,  F  (4  )  .  QLP(  20)  .  DT  ,  LF  ,HF  .  QP  ,  TS 

220 

COMMON  KPAR.KSTR.MVR.UMIN.UMAX.KPIT.KPIB 

221 

NTT = MTT- 1 

222 

101 

FORMAT ( 1 0F  8.3) 

223 

I F ( IKT-2  )  1.3,3 

224 

1 

REAO( 5 ,101 )XST 

225 

RE AD ( 5 . 101 ) YS 

226 

3 

J=1 

227 

DO  20  K=  . MTT 

228 

DO  10  N  =  J . 26 

229 

IF (XST(N)-XT(K) ) 1C.  10,2 

230 

10 

CONTINUE 

231 

2 

J  =  N-  1 

232 

L  =  N 

233 

YST(K )=( ( (XT(K)-XS  ( U ) )  T  (  Y  S  (  L ) - YS ( J ) )/ (XST (L)-XST (U) ) ) + YS ( J ) ) 

234 

20 

C  NT.-  :UE 

235 

YT(  1  )  =E(  1  )  'YST(  1  ) 

236 

DO  21  N= 1 .NTT 

237 

J  =  N+  1 

238 

21 

YT(J)=E(J)  (YST(J)-YT(N) )+YT(N) 

239 

RETURN 

240 

END 

END  OF 

FILE 

■ 


. 


...  QR<0,J<8  O  5 
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APPENDIX  F:  BOTTOM  PRODUCT  CONTROLLER  TUNING  BY  THE 
ZIEGLER -NICHOLS  CONTINUOUS  CYCLING  METHOD 
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The  distillation  column  was  operated  at  the  following 
condi t ions: 


Feed  Flow 
Reflux  Flow 
Steam  Flow 
Top  Composition 
Bottom  Compos i t i on 
Feed  Composition 


17.98 

g/s 

13.35 

g/s 

13.22 

g/s 

96.83 

% 

3.98 

% 

50.00 

% 

Various  proportional  gains  were  tested,  under  closed 
loop  control, for  a  2%  step  increase  in  set  point  to  obtain  a 
sustained  system  oscillation.  The  results  are  summar i zed  in 
Table  FI. 


Table  FI  Results  of  Bottom  Product  Composition 
Ziegler-Nichol s  Tuning  Runs 


Run 

Proportional  Gain 

System  Performance 

FI 

0.04 

offset 

F2 

0.06 

sustained  oscillation 

F  3 

0.12 

system  unstable 

F4 

0.20 

system  unstable 

The  control  and  manipulated  variable  responses  for  Run 
F2  are  plotted  in  Figures  FI  and  F2.  On  the  basis  of  the 
Ziegler-Nichols  recommended  settings,  the  controller 
constants  for  proportional  plus  integral  control  are: 

Kf  =.  0.45KPMy  =  0.027 

\  =  Pm/  1-2 


80.0  minutes 
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Set  Point  Using  A  Proportional  Controller  With  kd  =  0,06 
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Figure  F2  Steam  Flow  Rate  Response  For  Proportional  Controller  With  kp  =  0,06 

Of  Bottom  Product  Composition  Control  For  A  +2%  Step  Change  In 
Set  Point 
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APPENDIX  G:  SUMMARY  OF  PARAMETER  ESTIMATION  TESTS  FOR  STR 
CONTROL  OF  BOTTOM  PRODUCT  COMPOSITION 


The  closed  loop  parameter  estimation  runs  given  in  Table  G.1 
were  undertaken  to  identify  a  suitable  disturbance  flow 
pattern  for  obtaining  STR  parameter  estimates  for  use  in 
tests  of  the  regulator. 


Table  G.1  Summary  Of  Parameter  Estimation  Tests  For 

STR  Control  Of  Bottom  Product  Composition 

M  =  2  k  =  5  L  =  6 

^>c  =  10.0  P(0)  =  50.0 

For  stochastic  (Gaussian)  changes  in  feed  flow  rate: 


Run  Sland^nd.  InlLial 

Deviation  Eauamaieca 


EJLna.1 

Parameters 


JEanameter 

Behavior 


G1  1.0  g/s  0 . 0 


-0.2754  0.6673  Parameters 

-0.3796  0.1350  did  not 

-0.1461  converge 
0.0538 
-0.0746 
-0.0972 


G2  3.0  g/s  0.0 


-4.2410  -0.4138  Parameters 

-0.7210  0.6750  did  not 

-0.3219  converge 
0.5420 
-0.3189 
0.1649 


For  deterministic  step  changes 


21. 

3  g/s 

G3  10 

95 

,  Rn  700_ 

~ l _ 

^  nuu. 

14.4  g/s 


in  feed 

-1.2274 

0.2885 


flow  rate: 


0.2897  Parameters 

0 . 3 i 78  converged 

0.3419 

0.2607 

0.0360  * 

0.3481 


. 
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G4 


10 

18.0  g/s 


Final 

Parameters 


'14.4  g 


Of 

^00— lain.. 


Run  G3 


-5.5874  -0. 1607 
4.6377  -0.0247 
0.1027 
0.1419 
-0.0196 
-1.0294 


G5 


10  180 

IS 


18.0 


IP-Qmllb 

g/s 


14,4  g/s 


0.0 


-0.5038  0.0249 

-0.5259  0.2004 

0.3102 
0.2656 
-0.0625 
0.6574 


G6  10 


18.0  g/s 


700 
min . 


14.4  g/s 


Final 

Parameters 
Of  Run  G5 


-5.1508  -0.0412 
4.0556  0.0180 

0.0731 
0.0719 
-0.1195 
-0.9923 


d 


On  the  basis  of  these  results,  the  square 
isturbance  in  Run  G5  was  selected. 


Parameters 
did  not 
converge 


Parameters 

converged 


Parameters 

converged 


pulse 


' 


. 
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APPENDIX  H:  EFFECT  OF  SAMPLING  INTERVAL  ON  THE  STR  CONTROL 

OF  BOTTOM  PRODUCT  COMPOSITION 


The  number  of  parameters  to  be  updated  on-line  as  well  as 
the  amount  of  memory  required  with  a  real  time  process 
computer  are  directly  affected  by  the  k  parameter,  which  is 
defined  as  a  nonnegative  integer  expressing  the  system  time 
delay  as  a  multiple  of  the  sampling  interval.,.  Hence,  k  can 
be  manipulated  by  a  change  in  sampling  interval.  The  tests 
performed  using  different  sampling  intervals  are  summarized 
in  Table  H . 1 . 


Table  H.1  Summary  Of  FPC  And  PI  Tests  For  Bottom 

Product  Composition  Control  Using 
Different  Sampling  Intervals 
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By  increasing  the  sampling  interval  from  256  seconds  to 
512  seconds,  K  is  effectively  reduced  by  half.  The  FPC 
performance  is  superior  when  compared  with  the  PI 
performance  as  shown  in  Table  H.l.  When  the  sampling 
interval  is  further  increased  to  1024  seconds,  which  is 
equivalent  to  the  system  time  delay,  the  parameters  in  Run 
H12  were  still  drifting  at  the  end  of  the  run. 
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